This blank page was inserted to preserve pagination. 


SHAPE FROM SHADING: A METHOD FOR OBTAINING 


THE SHAPE OF A SMOOTH OPAQUE OBJECT FROM ONE VIEW 


Berthold K. P. Horn 


November 1970 


PROJECT MAC 


MASSACHUSETTS INSTITUTE OF TECHNOLOGY 


Cambridge Massachusetts 02139 


ACKNOWLEDGMENTS 


The author wishes to express his appreciation to everyone who 
contributed to the interesting research environment. 


‘Work reported herein was supported in part by Project MAC, an 
M.I.T. research program sponsored by the Advanced Research 
Projects Agency, Department of Defense, under Office of Naval 
Research Contract Number Nonr-4102(02). Reproduction in whole 
or in part is permitted for any purpose of the United States 
Government. 


This research was also supported by the South-African Chamber 
of Mines and the Council for Industrial and Scientific Research. 


SHAPE FROM SHADING: A METHOD FOR OBTAINING 
THE SHAPE OF A SMOOTH OPAQUE OBJECT FROM ONE VIEW* 


Abstract 


A method will be described for finding the shape of a smooth 
opaque object from a monocular image, given a knowledge of the 
surface photometry, the position of the light-source and certain 
auxiliary information to resolve ambiguities. This method is 
complementary to the use of stereoscopy which relies on matching 
up sharp detail and will fail on smooth objects. Until now the 
image processing of single views has been restricted to objects 
which can meaningfully be considered two-dimensional or bounded 
by plane surfaces. 


It is possible to derive a first-order non-linear partial differ- 
ential equation in two unknowns relating the intensity at the image 
points to the shape of the object. This equation can be solved by 
means of an equivalent set of five ordinary differential equations. 

A curve traced out by solving this set of equations for one set of 
starting values is called a characteristic strip. Starting one of 
these strips from each point on some initial curve will produce the 
whole solution surface. The initial curves can usually be construct- 
ed around so-called singular points. 


A number of applications of this method will be discussed including 
one to lunar topography and one to the scanning electron microscope. 
In both of these cases great simplifications occur in the equations. 
A note on polyhedra follows and a quantitative theory of facial make- 
up is touched upon. 


An implementation of some of these ideas on the PDP-6 computer with 
its attached image-dissector camera at the Artificial Intelligence 
Laboratory will be described, and also a nose-recognition program. 


4This report reproduces a thesis of the same title submitted to 
the Department of Electrical Engineering, Massachusetts Institute 
of Technology, in partial fulfillment of the requirements for the 
degree of Doctor of Philosophy, June 1970. 
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1, INTRODUCTION: 


1.1 SHADING AS A MONOCULAR DEPTH CUE: 


Consider a smooth object known to have a uniform surface. An 
image of such an object will exhibit shading (gradations of 
reflected light Intensity) which can be used to determine its 
shape, given only a picture froma single viewpoint. This is 
not obvious since at each potnt in the image we know only the 
reflectivity at the corresponding object point. For some 
points (called singular points here) the reflectivity does 
uniquely determine the local normal, but for almost alt 
points it does not. The shape of the surface cannot be found 


by local operations alone. 


Vy 
—& SOURCE 
NORMAL 


SENSOR 


Figure 2: Definition of the incident (i), emittance (e) 
and phase angle (g). 
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For many surfaces the fraction of the incident Light which is 
scattered ina given direction is a smooth function of the 
angles involved. It is convenient to think of the situation 
as depending on three angles: the incident angle (between 
local normal and incident ray), the emittance (or emergent) 
angle (between local normal and emitted ray) and the phase 


angle (between incident and emitted rays). 


It can be shown that the shape can be obtained from the 
shading if we know the reflectivity function and the position 
of the light-source(s). The reflectivity and the gradient of 
the surface can be related by a non-linear first-order 
partial differential equation in two unknowns. The recipe 
for solving this equation is to set up an equivalent set of 
five ordinary differential equations (three for the 
coordinates and two for the components of the gradient) and 
then to integrate these numerically along certain curved 
paths on the object called characteristics (5]. For while we 
cannot determine the gradient locally, we can, roughly 
speaking, determine its component in one special direction. 
Then taking a small step in this direction, we can repeat the 
process - the curve traced out on the object In this manner 
is called a characteristic, Its projection on the image 
plane will be referred to as the base characteristic. The 


shape of the visible surface of the object is thus given as a 
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sequence of coordinates on some such curves along its 


surface. 


An inittal known curve on the object is needed to start the 
solution. Such a curve can usually be constructed near the 
stngular points mentioned earlier using the known (loca! 
normal. The only addittonal Information needed is the 
distance to the singular point and whether the surface Is 
convex or concave w.r.t. the observer at this polint ~ such 
ambiguities arise in several other instances in the process 


of solution as will be seen. 


Figure 3: Image of a sphere and a stereo-pair of the 
characteristic curves obtained from the shading. 


To solve the equations, the reflectivity as a function of the 
three angles must be known, as well as the geometry relating 


light-source, object and observer. Multiple or extended 
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light-sources increase the complexity of the solution 
algorithm presented, But all of this initially needed 
information can be deduced from the image if a calibration 
object of known shape is present [In the same image. 
Furthermore, incorrect assumptions about the reflectivity 
function and the position of the ltight-source(s) can lead to 
Inconsistencies in the solution and it may be possible to 
utilize this information in the absence of a calibratton 


object. 


In practice it Its found that ff the object is at all complex, 
its image will be segmented by edges. Some of these are 
purely visual, due to the occlusion of one surface by 
another, others are angular edges (atso called joints here) 
on a single object. Another kind of edge Is the ambiguity 
edge. This Is an edge which the characteristics cannot 
cross, Indicating an ambiguity which cannot be resolved 
Locally. One can solve inside each region hounded by these 
various edges, but some global or external knowledge’ is 
needed to match up the regions. In the case of an angular 
edge on the object one can integrate up to the edge and then 
use the known location of the edge as an initial curve for 
another region (provided one resolves the ambiguity present 


here, as on all initial curves). 
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A very similar situation also obtains when one bridges a 
shadow. Since one edre of the shadow and tke position of 
the light source is known, we can trace atong the rays 
gBrazing the edge untii the corresponding image points fall on 
an illuminated region. Since we know the path of each ray, 
we can calculate the coordinates of the point where it 
impinges on the object. The edge of the shadow (which need 
not be on the same object) can now serve as an initial curve 


from which to continue the solution, 


A number of interesting applications of this method can he 
mentioned. The first of these concerns the scanning electron 
microscope (SEM) which produces [Images which are particularly 
easy to interpret, since the intensity recorded is a function 
of the slope of the object at that point and is thus a form 


of shading (as opposed to optical and transmission electron 


microscopes which produce intensities which depend on 
thickness and optical or electron density). The geometry of 
the scanning elect ron microscope allows several 


simplifications in the algorithm for determining shape from 
shading (e.g. there are no shadows), Because of the random 
access capability of the beam of this microscope fit should he 
easy and useful to combine it with a small computer to obtain 


three-dimensionat information directly. 
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Another important application ties tn the determination of 
Lunar topography. Here the special reflectivity functim of 
the material in the maria of the moon allows a very great 
simplification of the equations used in the  shape-from 
shading algorithm, The equations in fact reduce to one 
integral which has to be evaluated along each of a family of 
predetermined straight lines in the image, making for high 
accuracy. This problem was first tackled for areas near the 
terminator (the dividing tine between the illuminated and tke 
unilluminated part of the moon’s disk) by d. van Diggelen at 
the Astronomical Institute of the Netherlands [fn 1951 [2] and 
solved by T. Rindfleisch at the Jet Propulsion Laboratory in 
1966 [4] and the method applied to several alates returned 
by the Ranger spacecraft. This gave the first ftndication 


that the general solution discussed here might be possible. 


It should be pointed out that this method is complementary to 
the use of stereopsis, since the latter will match up sharp 
detail and edges while the shading information will determine 


the shape of the smooth portions of the surface. 


So far we have assumed that the surface is uniform In its 
photometric properties, Any non-uniformity will cause this 
algorithm to determine an incorrect shape. This is one of 


the uses of facial make-up; by darkening certain slopes they 
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can be made to appear steeper for example. In some cases 
surface-markings can be detected if they lead to 


discontinuities of the calculated shape. 


Judging by our wide use of monocular pictures (photographs or 
even paintings and woodcuts) of people and other smooth 
objects, humans are good at interpreting shading information. 
The short-comings of our method which are related to the 
shading information available can be expected to be found in 
human visual perception too. It will of course be difficult 
to decide whether the visual system actually determines the 
shape quantitatively or whether it uses the shading 
information in a very qualitative way only. A quantitative 
determination would involve operations more complicated than 
those used in edge-finding for example. Since the 
information is not local, the surface-shape calculations 


cannot be carried out entirely itn parallel. 


1.2 HISTORY OF THE PROBLEM: 


After formulating the image illumination equation as the 
basis of a method of finding shape from shading, a Literature 
search was performed to see if a solution had previously been 


obtained. The literature on perception has only ae few 
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conjectures on the possibility of determining shape from the 
monocular depth-cue of shading. Photogrammetry does not pay 
much attention to the reflectivity function, but only varlous 
integrals of it, measured by such devices as the integrating 
photometer. With few exceptions machine perception so far 
has been restricted to objects which can usefully be 
considered two-dimensional and objects bounded by planes 


(polyhedra). 


The one relevant research was found in the paper on tunar 
topography by T. Rindfleisch [4] which gives complete detalls 
of a solution obtained in the form of an integral in the 
special case of the reflectivity function of the moon. This 
raised the hope that a general solution existed. The 
(x’,y’,r) coordinate system used in [4] leads to Intractable 
equations - but we found aie solution ustng a different 
coordinate system, (x*,y’,z). As a check the solution for 
lunar topography was rederived from this set of equations 
(Rindfleisch found his solution in quite a different manner - 
searching for predetermined curves In the image along’ which 
the surface can be found as some integral tnvolving the 
measured image fllumination). A first program (old SHADE) 
was then written which solved along one characteristic at a 


time ustng various predictor-corrector-modifier methods [7]. 
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Another program (REFLEC) was used to measure the reflectivity 
function froma calibration sphere. Varfous short-comings of 
our image~dissector sensing device were affecting the 
accuracy of these measurements. Since very Little was known 
about the characteristics of this device on other’ than 
theoretical grounds [9], a program (TEXTUR) was developed to 
measure various propertles such as resolution, signal to 
nofse ratto, drift, settling time, scatter and pinholes in 
the photocathode. An attempt was then made to provide 
software to compensate for some defects such as distortion 
and non-uniform sensiftivity, using measurements from test 


patterns (DISTOR). 


These techniques allowed an estimation of what accuracy can 
be achieved under optimal conditions. The program had 
numerous problems when dealtng with objects other than sImple 
convex ones (mostly because it solved each charactertstic 
separately) and 6 must be apparent, was  sensftive. to 
Imperfections in the sensing device (partly because of the 


way {t obtained intensity gradients), 


After the defects In the first program had been found, and a 
decision made te rewrite it, a great simplificatton of the 
matn equatfons was found ustng a different coordinate system 


(x,y,z) anda slight extension of standard vector notation 
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(the voluminous equations for the inconvenfent coordinate 
system (x%,y%,z) are not reproduced here), An unfortunate 
but unimportant side-effect is an Increase in the comlexity 
of the derivation of the Lunar topography integral. The new 
equations and numerous changes in the method of solution were 
incorporated in anew program (new SHADE) which was less 
sensitive to the vartous shortcomings of our fmage-dissector, 
This program can handle objects somewhat more complicated 
than its predecessor and solves all charactertstics at the 


same time. 


In parallel with the programming work, theoretical efforts 
were made to define and get around some of the difficulties 
of the method of shape from shading. Of particular interest 
were applications where the equations simplify greatly. 
Unfortunately the massive simplification found in the case of 
lunar topography Is unique. Of most interest are the cases 
where we have some advance knowledge of the characteristics 
(for lunar topography they are completely Independent of the 
image - for the scanning electron microscope they are paths 


of steepest descent). 
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1.3 PREVIEW OF CHAPTERS - GUIDE TO THE HURRIED READER: 


References to articles and books listed at the end will be by 
numbers enctosed in brackets. Numbers contained in 
parentheses refer to sections and subsections in this work. 
In an attempt to be complete, a few subsections were Included 
which will have only timited appeal to some readers; hence 


this guide. 


Chapter 1 provides an introduction to the depth-cue of 
shading, its use in determining shape and its history. 
Chapter 2 develops the necessary equations in detail, 
Starting with the definition of the reflectivity function. 
Subsections 2.1.2 to 2.1.4 and 2.2 can well be skipped by the 
hurried reader, In section 2.3 the partial differential 
equation is obtained, the vector differentiation notation 
introduced and an equivalent set of five ordinary 
differential equations derived, Section 2.3 is perhaps the 
most important section. Sections 2,12 to 2.16 deal with some 
miscellaneous implications and may be omitted without loss of 


continuity. 


Chapter 3 describes in detail some practical situations where 
the special conditions encountered make use of the method of 


determining shape-f rom shading particularly attractive. 
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Section 3.7 deals with the scanning electron microscope. The 
reader should be warned about the tedious derivation of the 
simple integral for the case of tunar topography In 3.2 . 
Omitting subsection 3.2.3 will avoid the butk of the 
algebraic detail, and most of the conclusions will be found 


anyway in subsection 3.2.4 , 


Chapter 4 describes the experiments carried out with the two 
programs (using the results developed in chapter 2) to obtain 
shapes from images projected on an image dissector camera 
attached to the PDP~é computer in the Artificial Intelligence 
Laboratory. Section 4.1] deals with the less successful first 
program, and contains details on auxiliary routines. Section 
4.2 deals with the second program which solves the 
characteristics in parallel and also uses’ the important 
sharpening process, Sections 4.1 and 4.2 are next’ in 


importance to section 2.3 . 


Section 4.3 describes an application to a recognition task - 
that of nose-recognition., Section 4.4 contains an overall 
summary and conclusions about the capabilities of the method 
of shape-from-shading, with subsection 4.4,1 giving 
suggestions for future investigations. This is followed by a 


list of references. 
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2. THEORETICAL RESULTS: 


2.7 THE REFLECTIVITY FUNCTION: 


2.7.17 DEFINITION OF THE REFLECTIVITY FUNCTION: 


SOLID ANGLE dw 


SURFACE AREA dS 


Figure 4: Illustration of the variables used in the 
definition of the reflectivity function. 


Consider a surface element of stze dS inclined i w.r.t. the 
incident ray and ew.r.t. the emitted ray (The angles are 
measured w.r.t. the normal). Let the Incident Light 
intensity be I, per unit Seen perpendicular to the incident 
ray. The amount of Light falling on the surface element is 


then I, cos (i) dS. 


Let the emitted ray have intensity I, per unit solid angle 


per unit area perpendicular to the emitted ray, So the 
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vo RECEPTOR 


SVAFACE 


—- 
wees, 
~— 


cos(e)cos(g) - cos(i) 


tanta} Ss . COOeueesA gy ee 
cos(e)sin(g) 
cos(A) = cos(g) - cos(i)cos(e) 


sin(i)sin(e) 


Figure 5: Definition of the azimuth angle (A) and the 
projection of the emittance angle on the phase- 


angle plane -(a). 
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amount of Light Intercepted by an area subtending a_ solid 
angle dw at the surface element will be I, cos(e) dS dw. The 


reflectivity funcion b(i,e,2) is then defined to be I,/I,. 


If we want to be more precise about what units the intensity 
is measured in, we have to take into account the spectral 
distribution of the light emitted by the source, as well as 
the spectral sensitivity of the sensor (with this proviso we 
can speak of watts per unit area and watts per unit solid 
angle per unit area etc.). We need not be too concerned with 
this if we either use white paint, or measure the 
reflectiivity function with the same equipment later used _ in 
the shape-from=shading algorithm. It should be noted that 
for most surfaces the reflectivity function ts not 
independent of the color of the light used. Typically the 
specular component of the reflected light, being reflected 
before it has penetrated far into the surface, will be 
unchanged, while the matt component will be colored by 


pigments in the surface coating. 


Several other definitions of the reflectivity function are in 
use which are multiples of the one defined here by tr, 2, 
cos(e) and/or cos(i), The specific formulation chosen here 
makes the equation relating the incident tight intensity to 


the image illumination very simple. 
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2.1.2 FUNCTIONS DERIVED FROM THE REFLECTIVITY FUNCTION: 


The next few subsections (2.1.2.], 2.1.2.2 and 2.1.2.3) are 
Included to relate the reflectivity functions to those more 
commonly mentioned in the Literature. Some readers may want 


to skip these subsections. 


2.1.2.1 THE INTEGRATING PHOTOMETER: 
A flat sample of the surface under investigation is mounted 
in the center of a hollow sphere coated on the insitde with a 
highly reflective matt substance. Through one small hole a 
light ray enters and impinges on the sample with incident 
angle fi. A photosensitive device fs introduced through 
another small hole and measures an tntensity proportional to 


the light scattered by the sample into all directions. 


The total Intensity measured is: 
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SENSOR 


SOURCE 


A gonio-photometer (used for measuring 
reflectivity functions). 


Figure 6: 


INTEGRATING 
SPHERE 


Figure 7: An integrating photometer. 
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LT 11/2 
(ee dS sine) de dA 
O JO 
Vt Nar /2e 
eae [ff (i,e,8) },cos(e) dS sine) de dA 
O Oo 


Whe re cos(g) = cos(i) cosf(e) + sin€i) sine) cos (A) 


The total tncident intensity is I cos(i) dS, The fraction of 


light reflected is then: 
211 +T/2. 
bi) = f bCi,e,g)*(1/2)*xsin(Ze) de dA J/cos(i) 
ro) ° 


This function of the incident angle i has been measured for 
many paints and pigments, while the reflectivity function > 
is known for very few surfaces, Since it is difficult to 
relate measurements of I, to measurements of total reflected 
intensity, the device is usually calibrated with the sample 
replaced by a standard of known high reflectivity (e.g. MgO 
or BaSO, powder reflect more than 99% of the incident light 


in the visible spectrum). 
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Figure 8: Illustration showing quantities appearing in the 
integral for the integrating photometer. 


Figure 9: Illustration showing quantities appearing in the 
integral for the Bond albedo. 
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2.1.2.3 THE BOND ALBEDO; 
Another integral of the reflectivity function which is used 
is the Bond albedo, defined by astronomers as the ratio of 
total reflected light from a sphere divided by _the total 


incident light . 


If the incident intensity Is I,, then the ratio of reflected 


light to incident light is: 


tr/2 
C 1, [vo 2trr sin(t) cos(1) di J/(Iar’) 
oO 


Tr/2 
a b(t) sin(zi) di 
oO 
th SN 2 A W/2e 
. [ [ teivere sin(i) sin(Ze) de dA di 
oO oO (o) 


2.1.3 THE DISCRIMINANT 1+Z21EG-(I°+E 4G): 


Lom 
u 


In this subsection a discriminant Is developed which is 
needed in the program implementing the shape-from-shading 
algorithm. This section can be skipped without toss of 


continulty. 
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The three angles i, e and g, being the sides of a spherical 


triangle, have to satisfy the following relationships: 


ite 2g, etg >i and gti De 


It Is often convenient to express these three relattonships 


in terms of the cosines I, E and G of the three angtes. We 


first note that only one of the relationships could fail at a 


time. For example if ite < g : 


it2e < gte i.e. i < gte and 


2Zite < gti i.e. e € gti 


The angles are all posttive and tess than TT. Now assume 


that the condition ite < g holds, then: 


cos(Iite) > cos(z) 


since cosine is monotonic decreastIng for angles between 0 and 


TT. Expanding we get: 


cos(i) cose) - cos(g) > stn (i) sine) 


Since the right-hand side is positive, the left-hand side 


witl be too and we can square the expression. Using I, E and 
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G to stand for the cosines of i, e and g we get: 


CIE=<6)> > CPs 17 (71487) Tees 


1+21EG-(1° +E 74+G*) < 0 


We now have to prove the converse J.e. If the angles can 
indeed form a sphertcal triangle then the discriminant will 
be positive. Since i ¢ e+tg we have g ) i-e and similarly 
since e <¢ itg we have g }e-I, so: 

g > je-il and similarly | 3 [g-e| ande > |i-g| 
Applying the cosine as before we get: 

cos(g) - cos(1) cose) ¢ sin(i) sine) 
cos(i) - cos(e) cos(g) ¢ sin(e) sin(s) 
cos(e) - cos(g) cos(t) < sin(g) sin(i) 
From ite ¢ g etc. we get the same inequalities with the sign 
reversed on the teft-hand side. We needed to go to the 
trouble of showing that these inequalities hold for absolute 


values of the left-hand sitde, stnce fit no tonger Is 


constraint to be positive. So we have: 
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[IE-G] ¢ sin(t) sin(e) 


|EG-I| siIn(e) sin(g) 


€ 
|GI-E] ¢ sin(g) sini) 
Multiplying we get: 
| (IE-G)(EG-I)(GI-E)| ¢ (1-17)(1-E7)(1-G") 


Using one of the two signs for the right-hand side and 


expanding we get: 
(1-IEG) (14+21EG-(I7 +E74G*)) > 0 
Since [I], [E] and |G] ¢ 7 we have (1-IEG) > 0 and hence: 


1+2IEG-(174+E~4+G") > 0 


2.7.4 REFLECTIVITY FUNCTIONS AND THEIR MEASUREMENT: 


Surfaces where the three parameters i, e and g are not 
sufficient to fully determine the reflectivity are unsultable 
for this analysis (or at least reduce the possible accuracy). 
Examples are translucent objects and those with non-isotropic 


surface properties (e.g. the mineral commonty called tiger- 
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eye, hair, thin wax). Perhaps the most important determinant 
of the reflectivity functton fs the micro-structure of the 
surface (i.e. that structure smaller than the resolution of 
the sensor used in the determination of the reflectivity) and 
different reflectivity functions may apply at different 
magnifications (in addition, at high magnification objects 
become increasingly translucent). It ts best then to 
determine the reflectivity function under conditions similar 
to those used inthe determination of the shape of the 


object e 


One way to measure the reflectivity function ts to employ a 
gonio-photometer fitted with a small flat sample of the 
surface to be investigated. The device can be set for any 


combination of incident, emittance and phase angles, 


To avoid having to move the source and the sensor tnto all 
possible positions w.r.t. a flat sample of the surface when 
measuring the reflectivity function, it Is conventent to have 
a test-object which presents all possible values of i and e 
for a given g. (The constraints are ite ¢€ g, ete € i and 
gti ge). Use of such an object ts greatly simplified by 
using atelephoto lens anda distant source, giving almost 
constant g. It Is convenfent to tabulate the reflectivity 


versus f ande for each of a series of values of g. A sphere 
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is perhaps the eastest test~object to use tf one is willing 
to live with the decreasing accuracy in determining e as one 


approaches the edge. 


One could also have an object of known shape in the same 
scene as the object to be analysed. This solves the problem 
of having to know the source location and the transfer 
properties of the image forming system as well. In some 
cases objects of known shape and surface characteristics 
differing from those of the object under study are useful - 
for example a sphere with specular reflectivity can pinpoint 


the tocatton of the lLight~sources (e.g. the eyes in a face), 


2.1.5 MATHEMATICAL MODELS OF SURFACES: 


A number of attempts have been made to predict reflectivity 
functtons on a theorectical basis starting with some assumed 
micro-structure of the surface. White matt surfaces are 
usually finely divided grains of transparent material (e.g. 
snow and crushed glass), White paint usually consists of 
transparent ‘pigment’ particles (e.g. of S10, or TiO, ) of 
high refractive index (1.7 to 2.8) and small size (optimally 
about the wavelength of visible light) suspended in a 


transparent medium of low refractive index (1.3 say). If one 


PIGMENT TRANSPARE NT 
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Figure 10: Model of surface-structure. 


OPAQUE REFLECTIVE SUBSTANCE 


Figure 11: Another model used in the derivation of 
theoretical reflectivity functions. 
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choses a somewhat resulzr arrangement of suspended particles 
of uniform size and makes some very restrictive assumptions, 
one Con. derive: -e -petlectivity fenetier “and study- ats 
dependence on various parameters describing the model of the 


surface. 


Another type of surface is that of a highly reflective 
matertal (such as ametal) where the tight rays do not 
penetrate into the material. Choosing a particular type of 
surface depression anda statistical distribution of the size 
of these, one can again derive a reflectivity function. Only 
a few such models have been studied and little hope exists 
for modelling real surfaces well enoush and still deriving a 


closed expression for the reflectivity function. 


2.2 CALCULATION OF IMAGE ILLUMINATION: 


The equations derived in this section are only included for 
reference, since the program to be described later avoids 
their use by means of a normalization of the fmage intensity, 
These equations do have their importance however in 
justifying the choice of definition for the reflectivity 


function and in designing optical systems used in 


experimentation with the shape-from-shading method, 
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LENS IMAGE 


Figure 12: Diagram of optical system and quantities needed 
in the calculation of image illumination. 


Let d be the diameter of the pupil of the image forming 
device, f, its focal length and M the image magnificat ton 
(the ratio of the length of a tine in the image to the 


corresponding parallel tine on the object). 


Let a portion of the object surface of area dAg be Inclined 
at angle eto the tine from it to the imare-forming device. 


Its image will have an area of dA, = dA, M /coste). 


Let the incident intensity at the object patch be I, per unit 
area perpendicular to the incident ray. Then the emergent 
intensity per unit solid angle will be I, = 1, $(1,e,%). The 
light captured by the image forming device is I, dA dw/cos(e) 


where dw is the solid angle formed by the cone of angle a. 
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dw = 211 (l-cos(a/2)) = 41t1sin (a/4) 


We would Like to express this in terms of M and the so-called 


f-number fr: 


fn = 1/(2 sin(b/2)) 


Ji / 4) (14h)* (Fo /d)* = (14M) (fe/d) if (fo/d)>1 


The f~number usually indicated on the lens-barrel is (f,/d) 


or /(1/4)+(fo/d)* . 


fFr-(1/4) = (14M)? (f,/d)™ 
cos(@/Z) = eS 
/(M 7/4) +(14M)* (Fo /d) 
(4fe-1) 
cos (Q@/2) = SSS 
(4f,4+M*-1) 
Me 
dw = 217 (1-/1- —————_ ) 
(4fR+M*-7) 
Mz 
> ta —————_ if f,x> 7 


(4f24M*-1) 


The intensity per unit area in the image is: 
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I, [dA,/(cos(e) M*)Jxfcos(e)/dA ] dw 


— 
wy 
t 


I, dw/M* = I, 2n/(4f,) 1f M <1 


R 


I, = 1,6 dw/M* I, @n/C4fp) tf M <7 

It becomes apparent why we chose to define the reflectivity 
function the way we did and also why one might want factors 
of tr and/or 2 in the definition. It should be noted that In 
practice one does not usually employ this equation, but 


rather normalizes the expressions used, 


2.3 THE IMAGE ILLUMINATION EQUATION: 


This section contains the derivation of the Image 
illumination equatton and the analytical formulation of the 


shape-from-shading problem. 


2.3.1] PREVIEW OF HOW TO OBTAIN THE PARTIAL DIFFERENT IAL 
EQUATION: 


At a known point on the object we can calculate g. We would 
Like to find the gradient (or at least its component in one 
direction) at this point so as to be able to continue the 


solution to a neighboring point, Measurement of the Light 


SOURCE 


Figure 13: 
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yas 


Details of the geometry of image i]lumination 


and projection in the imaging system. 
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reflected tells us something about | and e. Since only one 
measurement ts involved, we cannot generally hope to 
determine both i and e tocally, but only a relation between 
the two. There are exceptional points where the normal is 
locally fully determined and this ts useful in finding 
initial conditions as explained later. This situation ts 
contrary to that obtaining in the use of texture gradients 
(see section 2.176) where the gradient fs known locally 
(except fora tworway ambiguity). In obtaining a solution 


from the shading, only a global approach will work, 


Collapse the two principal planes of the image-forming system 
together, forming the x-y plane. Let the z-axis coincide 
with the optical axis and extend toward the object. Let f be 
the exit pupil to image plane distance and assume that the 


image and object space refractive indexes are equal. 


Let t be the ratio of image illumination to object luminance 
(can be found from taws of optics - see section 2.2). Let 
a(x,y,z) be the Incident Light intensity (usually constant or 


obeys some inverse square law), Let A(x,y,z) = t*a(x,y,z) 


Let r= (x,y,z) be a point on the object and is = OK GF) 


the corresponding point in the Image. 
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b(x’,y*) is the intensity measured at the imasze point 


(x*,y”). 
Let I = cos(i), E = cos(e) and G = cos({ge) . 
We have A(r) $(1,€,G) = b(r*) 

w ww 


Let p anda be the partial derivatives of zwer.t. x andy. 
We would like to show that this equation involves x, y, Z, P 


and q only. 


2.3.2 NOTATION FOR VECTOR DIFFERENTIATION: 


If A is a vector (3-tuple), then A = [A| is the magnitude of 
A. Also tet AK = A/ A be the corresponding unit vector. 
Consider the dot-product A.B as matrix multiplication of the 
1 by 3 matrix A by the 3 by 7 matrix BY (the transpose of B), 
Consider partial differentiation w.r.t. a vector as the 3- 
tuple whose components are found by differentiating w.er.t,. 


each component fn turn. Then for example: 


> 
u 
>> 
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At times we will also need the parttal derivatives of vectors 
Wer.t. vectors, These are defined as 3 by 3 matrices (the 
first row being the result of differentiating w.er.t,. the 


first component and so on), then for example: 


We will also use parttal derivatives of dot-products of unit 


vectors w.r.t. vectors. For example: 
an 
X = A.B and we want Xp, 
e e s “ “a 
To avoid finding A, we write A X = A.B and then: 
A,X + A Xa = A 


e 
& 


Extending the definition of dot-product in the appropriate 


way we find: 


in PD ON xe " 
A .B={010]B = B 
~S 001 ss 
AX, =B-xXA 


~< 
> 
tt 
-o 
~ 
> 
pa 
2m> 
t 
>< 
Ur) 
eet 
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2.3.3 THE EQUATION IS A FIRST-ORDER NON-LINFAR P.D.E.: 


If the reflectivity function is b(1,£,G), the normalized 
incident light intensity at the point r = (x,y,2) is A(r) and 
the intensity at the corresponding Image potnt a = (x’,y’,f) 


is BOE) then: 
A(r) O(1,E,G) = b(r*) 


This ftmage fllumination equatton fs the main equation studied 
here. When finding a solution we assume A(r) and $(1,6,G) 
are known and bCr’) ts obtained from the image. We want to 
show that the equation is a first-order non-linear partial 
differential equation in two Independent varfables, i.e. an 


equation of the form: 
F(x, y, #, P, QQ) = 0 
where p = z, and q = zy are the partial derivatives of 2z 


weret. % andy respectively. From the simple projection 


geometry we have: 
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where f is the exit pupil to imare plane distance, We took 
. « . ¢ eo 

care of image reversal by orienting the x and y axes 

appropriately. It remains to show that 1, E and Gare 


functions of x, y, #, p and q. An inward normal to the 


surface at the point EAS Ae ip gy. 1.3 
Let the light-source be at ae (Ke Ve PEs 21 Then the 
incident ray will be Or ae ees ee ar and the emergent 
ray “Pe F ovr. Clearly then: 
nw “aA “aA an n~ “a 
TS eg: EO oa re San Ge ee Rae 
nN 4 . 

Where the s denote unit vectors. All the terms thus 


involve only x, y, #, p and q. It follows that we are dealing 
with a first-order non-linear partial differential equation 


in the two unknowns x and y. 


2.3.4 SOME DERIVATIVES NEEDED IN THE SOLUTION: 


When solving the P.D.E. by the method of character!stics we 
will nmeed the following partial derivatives (see section 
2.5), which it is convenient to introduce here, following the 
expansion of I, E and G in terms of dot~products. Using the 


results developed in subsection 2.3.2 we get: 
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described earlier (values of x, y and z) plus the values of p 


and q on them, 


Since we can multiply the eauation F = 0 by any non-zero 
smooth function A(x, y, Zz, Pp, a) without altering the 


solution surface, we can obtain a different set of equations: 


Ke Neg y =), z= (pF, + af, ) 


p = \ (- Fy "7 pF; ) and q = Xe Fy a aF, ) 


The solution to this new set of equations will differ only In 
the values of the parameter $s at any given” point. For 


example if we let 


Poa, 5 


= 11 /F% + iy SCOR gb 


the parameter Ss gives us arce- length along the 
characteristics. This is used In the programs to. be 
described later. Of course we can only do this if the 


denominator is not zero; at singular points and ambiguity 
edges It will be zero (t.e. Fp = Fa = @ and A%00), A 
different choice for ® will be used tater in the discusston 


of the scanning electron microscope (section 3,1). 
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2.3.6 OUTLINE OF PROOF OF EQUIVALENCE OF THE SET OF 


O.D.—.°S TO THE P.D.E.: 


In this subsection the equivalence of the five ordinary 
differential equations to the image illumination equation is 
discussed. The reader who belteves the equivalence hotds may 


well skip this subsect ion! 


At a given point (x9, Yo, #,) the equation F(x, y, 2, Pp, q) = 
0 represents a relutfon between p and q. That fs, it confines 
the possible solution normals at this point to ae one- 


parameter family of directions [5]: 
(-p(%), q(m%), -1) 
Increments in the feasible tangent planes thus satisfy: 
dz = po(x) dx + q,(%) dy 
Differentiating w.r.t, &% we get: 
0 = pila) dx + ql(m) dy 
(Dashes are derivatives w.r.t. &% In this subsection). But 


by differentiating the equatifon F(x, y, z, p, q) = 0Ow.r.t.& 


we also get: 


Page 47 
F,pa(®) + Feast ) = 0 


Hence: dx/Fp = dy/Faq 


What we need now are stmilar equations for feasible 


increments in z, p and q. First we have: 

dz = p dx + q dy 
in the solution surface (this surface is selected from all 
possible ones by choosing one passing through a given Initial 


curve - see later, subsection 2.3.7). Hence: 


dz/(pF, + qF,) = (p dx + q dy)/(pF, + aFq) 
= dx/F, = dy/Faq 


Finally differentiating F(x, y, #, Pp, q) Ow.r.t. x and y: 
Fy + Fup + Fon. + Faq, = 0 
Fy + Fea + Fep, + Faay = 0 

dz Oz 


but Py = Gy = 


Oxdy Dydx 


dp/ (F, +p Fa ) ~ (p,dx + py dy)/(p,F, + pF, ) - dx/Fp 


daq/(F, +qF) 


u 


~ Ca,dx + aydy)/(a,F, + ayF,) = ~ dy/Fa 
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dx/Fp = dy/Fa = dz/(pFp + qF,) 


~dp/(F, + pFa) = - da/(Fy + aF3) 


Introducing the parameter s we get the five O,D.E.’%s 
mentioned earlier. We have shown that a solution to the 
P.D.E. must also satisfy theses five 0.D.E.%s, It is a bit 
more difficult to show that a solutton to these 0.D.E.’s ts 
necessarily a solution of the P.D.E. (see [5], page 28). 
Basically it needs to be shown that the equations for p and q 
produce results which continue to be consistent (l.e. equal 


to the partial derivatives of z w.r.t. x and y). 


2.3.7 INITIAL CONDITIONS NEEDED: 


To select a particular solutfon surface amongst all posstble 


solution surfaces one needs to specify an fnitial curve 


through which the solutfon surface must pass: 


x = x(t), y = y(t) and z = z(t) 


Along this curve we must satisfy: 


INITIAL | 


CURVE >? we 


Figure 14: Characteristic strips sprouting from an 


initial curve. 


Figure 15: Illustration showing quantities appearing in the 
integral for the case of uniform illumination 


(Similar to Figure 8). 
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= 0 and for a truly distant source: 


Replace an by kr. and let k-»00 then: 


I, = 0, 1, unchanged 


E,= 0, E, unchanged 


ep) 
u 


~ * Org ME; - GT), Gy * 0 


In addition choosing the z-axis along rm; removes further 
terms. 
2. DISTANT CAMERA: (Telephoto lens or the object subtends a 
small angle at the camera) 


Replace re by kr, and tet k->oa then: 


I. and I, unchanged 


FE, = 0, E, unchanged 


<p) 
1 
ul 


C/red (Fe - GT), Ga = 0 


o~ 


In addition choostng the z-axis along re removes further 


terms. 
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3. DISTANT SOURCE AND DISTANT CAMERA: 


I~= 0, I, unchanged 
Er 0, E, unchanged 
Gr = 0 , Gy = 0 


Most practical situations are an approximation of this 


case. 


4, SOURCE AT THE CAMERA: 


he. ue I= —E andG =I] 
I, = Ey, unchanged 


I, = Ex unchanged 


> 
T 
u 


0 and Gy = 0 


5. DISTANT SOURCE AT DISTANT CAMERA: 


E, = G, = 0 


r 


= 
a+ 


_ 
f 


ws 


= E, unchanged, G,= 0 


Choosing the object to be on the z-axis removes’ further 


terms. This is the simplest possible case. 
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6. UNIFORM ILLUMINATION: 


Uniform itluminatton (or an approximation thereof) Its 
fairly conmon and might at first sight appear not to fit 
into our framework, This subsection shows the 
equivalence of uniform flluminat fon to one where a point- 
source is at the camera anda different reflectivity 


function obtains. 


The integrals here are analogous to the ones obtained for 
the integrating photometer except that we have constant 
emittance angle rather than constant Incident angle. If 
the Incident light intensity Is I, per unit area oriented 
in any directton, then it {fs easy to show that I,/m falls 
per unit solid angle per unit area perpendicular to It. 
The emitted light per unit solid angle per unit area 


perpendicular to the emitted ray fs thus: 


Ie Atr/2t 
I, Ve) = ra reoef [gttrererecrreresinct di dA J 
ro) °o 


/cos(e) 


This is the same stItuation as ff we had a source at the 


camera anda reflectivity function such that: 
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(Except that for uniform {tllumination a certain amount of 
self-shadowing can occur for non-convex objects) 
2.5 THE FIVE O.D.E.°S FOR THE IMAGE ILLUMINATICN EQUATION: 
F(x, ¥, @ P, a) = ACr) $CI,E,G) - b(r*) = 0 
We know A(r) and b(1,E,G), and obtain b(r’) from the image. 


We need Fy, Fy, Fa, Fp and Fa. Since r = (x,y,z) and 


ns (-p,-q, 1) we can get all of these derivatives from Fr 


and F,. 
Fe = ACr) O.C1,E,G) + Ay Cr) @CILE,G) ~ by (17) 
Fr = ACr) $ (1,E,6) 

Let = (1,E,G) then: 


a 
dr $02 


It 
a] 
3 
a 
~o— 
2s 
fl 
Pee = oe 
20 
det] 
es 


Note that a, and a, are 3 by 3 matrices, the rows of which we 


~ 


computed in a previous subsection (2.3.4). 


rae) 
24 


Note that 


Similarly: 


a.) 
Is 


FI G=A1 


c 
al 


= | (i/rea@- 


WW Ia PD 
= 
st 


e 


) + In Me-6 B,) 


1/r; -I/ry 0 n 
= 1/¥ 0 ELCs e. 
0 (M/re-G/re) i /r,-G/rg)]} \ EF, 


this is the product of two 3 by 3 matrices. 


(in) (re - Tn) 
Ci/n) CF, - & 


0 
-I/n 1/n 0 n 
= | -E/n 0 1/n fr; 
0 0 0 (a 
ne 


The calculation of b (r’) from br”) will be described in 


section 2,7 
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2.6 CAMERA PROJECTION FQUAT IONS: 


The projection equations derived here are used in section 


247 


So far we have assumed the camera to be at the origin 
oriented with its optical axis directed along the z-axis and 
the image-plane x* and y” axes parallel to the x and y axes, 
Moving the camera fromthe origin tntroduces only a minor 
change in the equations, If however the camera is oriented 
ina different way, some of the equations become more 


complicated. 


Let R be the orthonormal 3 by 3 matrix which takes the z-axis 
into the optical axis and the x and y axes into the x”* and y” 


axes. Then: 


whe re 2 = (0,0,1) is the unit vector along the z~axis In the 
image coordinate system. ir is the vector fromthe exit- 


pupil to the image point in the same coordinate system, 


If two images are taken with the camera orlented differently, 


the area recorded in both Images will be spatially distorted 
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only. That is, a simple transformation will take the one 


image into the other. 


Where R, and R, are the two rotatfon matrices and f, and f, 
the corresponding exit pupil to Image plane distances. This 
transformation is useful if one wishes to orlent the optical 
axis alone xr, or r, (to simplify the equations for the 
derivatives). 


2.7 OBTAINING INTENSITY GRADIENTS: 


To evaluate the derivative F, (section 2.5) we need b,(r’). 


BG Dis Di dace: sh 
f f Rr oz 
Pe sdk) iene = a CR) 
ww ~ (Rr). (Rr.z)* 
f (Rr) (RZ) 
mo (R -  ——————_ 
(Rr) .2 (Rr) .2 


In the sItmple case that the camera is orlented properly: 
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om SDS 


Sal 


1 
= =-1{ 0 
z\0 


Written out In full we have: 


an OS 
ae | 
<x 
on“ ™ 
TH IW 
we 


(b,, by, b,) = (f/z)[b,,, by, ~((x/z)b,. + (y/z)by-)) 
by. and by- are measured directly from the image. 


Since the intensities measured from the Image do not locally 
determine the normal, one might well ask what, roughly, such 
measurements do determine. The components of the gradient of 
the Intensity are related to the second derivatives of the 
distance to the surface, while the Intensity itself is 
related to the magnitude of the first derivatives. This 
relationship becomes exact for the case of a distant source 
at a distant camera (section 2.5, case 53; see also 3.1.2). 


It should be noted that the equation for F, (section 2.5) 
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also Involves A,. Usuatly A is fairly constant over the area 
of the object recorded in the image, or at least satisfies a 


simple inverse-square equation, 
If A= (r,/r-), then Ay = -2Cre/rt ry 
= r re ? en r a re re re . 


Where rj, is the incident vector, and re Is the length of the 


incident vector to the singular point. 


2.8 OBTAINING INITIAL CONDITIONS: 


It would be a great disadvantage if one always required = an 
initial curve to start the solutfon from. Fortunately it Is 
usually possible to calculate some initial curve if one makes 
some assumptions about the surface and uses the. special 
points where the reflectivity uniquely determines the local 


normal - these points will be called singular points. 


2.8.1 USE OF THE SINGULAR POINTS: 


The singular points are the brightest or the darkest points 
(depending on the reflectivity function). At all other 


points the normal cannot be locally determined. The singular 
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points are points corresponding to values of i ande_ for 
which the reflectivity ts a unique global maximum or minimum, 
These may be either extrema in the calculus sense or at the 


Limiting values of the angles. 


This method cannot be used ff the surface does not contain a 
surface element oriented in this special direction. The 
points are found by looking for the brightest (or darkest) 


points in the image. 


Alt we still need to know then fs the distance of this point 
from the camera, but since one is usually only interested in 
relative distances this Is mot a_ serious restriction. 
Unfortunately it wltl be found that the solution will not 
move from these singular points, le. x*° = y* #0, This Is 
an indication that the algorithm needs to be tnformed about 
which way the surface is curved (convex or concave). To make 


this more concrete assume we have a_ distant source and can 


thus calculate G at each image point. 
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2.86.2 THE SOLUTION WILL NOT MOVE FROM A SINGULAR POINT: 
Consider the variation of > with — firsts: 
1. If the extremum occurs for O0<E<? then b. = 0, 


2. If the extremum occurs for E = 1 then f = f, and hence 


E, = (nrg - EA) = 0. 


3. If the extremum occurs for E = 0 then ne = 0 and 
En = (/n)k. That Is, xpt+yp-z=0 and Ep = (1/nr)x and 
Eq coal (l/nrdy. 
x = o-(1/nr)x and y = de(1/nrdy 
zB = oe(1/nrd(pxtay) = o.(I/nrdz 


In case I and 2 we have by E, and $e Ey = 0, 


Now constder the vartation of $ with I: 


ul 
S&S 
. 


1. If the extremum occurs for 0<I<7 then >. 


hence 


3 
tt 
Yagy 
fel) 
3 
om 


2. If the extremum occurs for I = 71 then 


T,= (/n(e, - TT) = 0. 


3. If the extremum occurs for I= 0 then 
I,= (1/n)ty. That is (x-x,)p + (y-y,)q ~ (2z-2,) = 0 and 


Ip = (I/nre)(x-x5) and Iq = (I/nrj)ly-yo). 


Me 
i] 


b (1/nr,) (x-x,), y = ¢ (l/nr;) (y-y,) 
@ (1/nr,)((x-x Pp + Cynye)a) 
> (I/nr, )(a-2,) 


Ne 
i] 


In case | and 2 we have $s I, and os I, = 0, 
Now x = Fp and y = Fy, and zs pF, +taFa. 
Ar) $(1,E,6) 


$y], a >. Ep 
$s = o,1, + PE, 


mn 
3 
" 


o- 
i 
" 


So in all combinations of cases I and 2 for FE and cases 7| and 


2 for I we find x = y = 0 and hence also # = 0, therefore: 


That fs, the projection of the solution point Into the image 


Is not moving as the parameter s fs changed. 
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In the case E = 0 we find that even though x and y may be 
non-zero, x° and y” = 0, 
All that remains is the case I = 0. Here too x’ = y’ = 0, 


when the source is at the camera or if — is undetermined, We 
have found no paints with an extremum for I = 0 where E was 


determined (i.e. the global extremum was not unique), 


2.8.3 GETTING THE INITIAL CURVE FROM A SINGULAR POINT: 


If the surface is convex (or concave) at the singular point 
and we have a guess at the radius of curvature (from the 
overall size of the object for example) we can get around the 
problem of singular points by constructing small spherical 
caps on them. Difficulties will be encountered if this point 
happens to be a saddle point (The presence of a saddle poltnt 
however usually Indicates that other singular points exist 


where the surface is either convex or concave). 


Let s be the vector from the camera to the singular point 
(found from its known image coordinates and its distance from 
the camera). Ris the estimated radius of curvature and / 
the distance we decide to step away from tte singular point 


(determined in practice by considerations of uncertainty In 
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Figure 16: Construction of the initial curve near a 
singular point. 


SINGULAR 
POINT 


Figure 17: Illustration portraying three solutions obtained 
for varying initial radius of curvature - 
showing the small effect which errors in the 
initial curve have on the solution. 
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the position of the stngular point and the desired detail in 
oa 
the solution). The known normal at the stngular pofnt Its Noe 


“Aw 
We construct a spherical cap with center S - RN. 


Let RV = RY - ¢ 


A 
Si = $+ (RRM 


Let X= ax fh, where y = (0,1,0) 
Y= fixX 
T(t) = p IX cos(2mt) + y sin(27r t)] 0gt<l 


Points on the initial circle are then given by 


$, + T(t) 
We also need an initial guess at p and q, so we construct Nv, 


(an outward normal): 


The requirement for an initial guess at the radius’ of 
curvature jis not as restrictive as It might seem, since the 
required accuracy is extremely low, This ts because is 
usually very much smaller than R, and hence a_ change in R 
affects the position of the inittal curve very little. Even 


more tmportantly the values derived for p and q need not be 


accurate since they are only used as a first guess in an 
iterative method of finding p andaq on the initial curve 


before starting the solutton. 


2.9 NON|=POINT SOURCES: 


Uniform scurces have already been dealt with. Perhaps the 
easiest other case is acfrcularly symmetric source at a 


distance large compared to the dimensions of the object. 


269.4 CURC ULARLY SYMMETRIC “SOURCES4 


Distant circularly symmetric sources can be replaced by a 
point source after modifying the reflectivity function. One 
merely convolves the reflectivity function with the spread 
function of the source (a bit of spherical trigonometry is 
involved here). Strictly speaking one should perform the 
same operation with the entrance pupil of the camera since it 
too subtends a finite angle at the object and accepts a 
bundle of light-rays. Since > Is smooth (except at I = 0 and 
I = 1) it will be changed very little except at these points. 
The main change will be that ¢ does not tend to 0 as I. tends 


to 0, but rather for some negative value of I. Also the 
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specular component will be more smeared out. 


CENTER OF 
SOURCE 


OBJECT 


Figure 18: Illustration of circularly symmetric source and 
quantities used in the convolution. 


Let the source Intensity be I(a) per unit solid angle at the 
angle a from its center when viewed from the object. Then 


the new reflectivity function $°(1,£,G) is? 


21 A Ao 27, Ao 
four, I(a) $€1°,E£,G*) a da dv f fics a da dv 
° ° . ° 9 


Where a, is the total angular diameter of the source, 
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And cos(A) 


(cos(g)-cos(i) cos(e))/(stin¢(i) sin le)) 


cos(i*) = cos(i) cos(a) + sin(i) sin (a) cos(v) 
sin(SA) = sin(i*) sin(€a)/sin(v) 
cos(s”) = cos(AtSA) sin(i*) sinle) + cos(1%) cos (e) 


2.9.2 MULTIPLE SOURCES: 


When the source distribution is not easily treated as above 
one can introduce a different Ay for each source and replace 


the main equation by: 
DAgGE) PCy LEG) = bE) 
k 


Difficulties in finding init fal condit fons witt be 
encountered with multiple sources unless they are of special 


kinds (e.g. a point source and a uniform source). 


2.10 TYPES OF EDGES: 


Several kinds of edges appear in an t!mage - each with fits own 


properties and problems for our atgorithm: 
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Overlap - (occlusion of one object by another) 
discontinuity in z. The program must detect this or It 


witl erroneously continue a solution across such an edge. 


Joints - (angutar edges on an object) discontinulties in 
the derivatives of z. One cannot continue p and q_ across 
such an edge. It Is possible however to use the position 
of the edge as a new Initial curve. This and the previous 
condition can be detected as a step in the fntensity 


distributfon or froma highlight on the edge. 


View edges - special case of 1. , where no joint appears, 
i.e. the surface is smooth and E tends to 0 as we approach 
its This is easily detected by the program during the 


calculation of the solution. 


Shadow edges - here I tends to 0 as we approach the edge 


and again the program can easily detect this. 


Other edge of shadow - if the shadow was bridged this edge 


may serve as a new initial curve. 


Ambiguity edges - some are tines of aggregation of 
singular points (on which »~7 ys The characteristics 


will not cross an ambiguity edge (see sectton 3.1.3). 
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2.717 SHADOWS AND SEL F-ILLUMINAT ION: 


If the single source is not at the camera, shadows will 
appear. Soluttons can be carried across shadows since the 
position of the source is known and one can construct a ray 


through the last [Tlluminated point and trace it until it 


meets another illuminated region. Only the coordinates and 
not the local gradient of this new point will be known, It 
is necessary to carry this operation out for all 


characteristics entering the shadow, producing a new initial 
curve at the other edge of the shadow where we can restart 
the solution. In practice care has to be taken because of 


noitsyness of the solution. 


Self-illumination is a difficult problem to deal with unless 
the object is convex or [ts albedo is low (less than 0.4), 
An estimate of the effect of self-illumination can be 
obtained from a consideration of two semi-infinite matt 
planes joined at right angles. These are illuminated from a 
very great distance and such that the incident rays make an 
angle a with one of the planes. Let the reflectivity of the 
surface obey lamberts law and the fraction of the Incident 
light reflected be k. Contrast between two Intensities I, and 


I, is usually defined to be: 
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Figure 19: Bridging a shadow. 


Figure 20: Two semi-infinite planes joined at right angles. 
Used in the study of self-illumination. 
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If we ignore light reflected more than once, we find the 


contrast between the two planes to be: 


C, = tan (a - 11/4) 


White if the self-illuminatton is taken Into account we get: 


2 =k 
C, = tan(a ~ tr/4) 
2 +k 


Contrast is thus reduced by a factor of (2 - k)/(2 + k). 


This factor vartes from 1/3 to 7? as k varles from 7 to 0, 


Note: the rest of chapter 2 contains some miscellaneous 


items that did not fit In elsewhere, 


2.712 THE INVERSE PROBLEM - GENERATING HALF-TONE IMAGES: 


The inverse problem of producing images of a specified scene 
with shading and shadows Is vastly different from the method 
of shape-fromshading. Most programs written for this 
purpose can be used for objects bounded by planes only. The 


main tssues of optimization of the calculation of which 
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surfaces are visible to the source and camera respectively 
have been dealt with in some detail in recent work [8]. 
Although the two problems are inverses of one another, the 


methods used are quite different, 


An interesting problem of a mathematical nature (and 
Incidentally with application to cutting wood-cuts) is that 
of producing curved Lines In a plane such that the density of 
lines is proportional to the shading in the image of some 
real or imagined object. Preferrably one would tike as small 
a number of “unneccesary” breaks in the Lines as possible, 
i.e. the lines should either close on themselves or leave the 
image. Another restriction one might apply is that the lines 
should not cross (When producing wood-cuts one would most 
likely also reflect some of the surface texture in the cholce 


of lines). 


For a special case, a solution is immediately at hand. This 
ts the case where we have a distant camera at a distant 
source (section -2.4, case 5; see also 3.1.2) and a 


reflectivity function b such that: 


OC1,1,1) = 1 = 1/ /14+p*+qt 
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Here the contour lines give a solution, with no crossing 
lines and no “unneccesary” breaks, One of the most 
attractive feature of contour maps is perhaps just this fact 


that they provide some shading information. 


2.13 HUMAN PERFORMANCE WITH MONOCULAR PICTURES: 


Jugding by the popularity of monocular pictures of people and 
other smooth objects, humans are good at interpreting shading 
information, Since they use the same bastc information as 
our shape-from-shading algorithm we expect to find similar 
short-comings (see section on facial make-up for examle). 
Supposing the human visual system does not use the’ shading 
information in simple heuristic ways only, one might expect 
that the perception system “solves*® the equations or a much 
simplified form of them. Since this cannot be done tocally 
(the way some portions of an edge-finding process might work) 
it ts difficult to suggest an elegant and sImple mechanism 
and a place to look for it. Presumably ft would have to 
Involve computational waves travelling outward from the 


singular points. 
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2.14 ERRORS AND INCONSISTENCIES: 


It ts difficult to estimate analytically the error In the 
solution because the equations are so non-linear. d, b, and 
A cannot be measured to better than 5 or 10% accuracy and 
numerous practical problems such as non-uniform sensitivity 


of the sensor have to be taken care of, 


Only a simple error analysis can be presented here. Suppose 
we wish to determine the effect of varying Inclinations on 
how a given error in the input data (Intensity In the image) 
relates to errors in the coordinates determined on the 
characteristics. We need to determine the rate of change of 
p w.r.t. b. Consider a particularly simple case, that of a 
distant source at a distant camera (As has been’ mentioned 
previously and will be demonstrated in section 3.3, the 
equations for this case are particularly simple). Next 


assume that one of the gradlent components, q say, Is 0. 


We have b/A =.O(1) = $(1/ /T+p24q> 2 = $C /14p™ ) 
Then p = /1/(o7' (b/A))* - 1 


We need to differentiate p w.r.t. the ratio lL = b/A. 
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; -b7 (1) 


Pe ees | hss eee 
t CL) /1= 0b 


For both I-— 0 and I — 1], the error In p becomes very large 
for a glven error in t (since in the first case q-"C1) —» 0 
and in the second case o-'C1) ooo ss a a This ts not very 
surprising since in the first case we are looking 
perpendicularly down on the surface and I will vary very 
slowly with p, while in the second case we have near 
tangential incidence and smalt changes in the angle of 
Incidence (and hence also I) will correspond to large changes 


in p. 


We note that In this rather special case, the error 
contribution to the solution ts large in some areas, while 
being smalt in others where the incident angle is not to 
close to 0 or T/2, The actual error will also depend on qo! 
and the error in measuring b/A. In a case with less 
restricted tighting condittons the relationship between the 
inclination of the surface and the error-rate will be more 


complex. 


We considered the derivative of p w.r.t. lt, since ft fs the 
integral of the error in p which constitutes the error [In z 


for any one charactertstic. 
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e(s) = Sp(s) ds = p §1(s) ds 


Where e(s) is the error in z for a given characteristic as a 
function of are-distance from the singular point, p(s) Is 


the error in p and {l(s) ts the error in Ll. 


In this context one may also want to discuss inconsistencies 
in the solution. If either the Lighting conditions or the 
reflectivity function are incorrectly specified, an Incorrect 
shape will be calculated. The shape determined may or may 
not violate the requirement of smoothness. If the calculated 
shape Is not smooth it can be concluded that the solution (at 
least [In some region) is incorrect, and that the given source 
positton or the given reflectivity function are incorrect. 
It is easy to give examples of the case where false 
assumption will lead to a smooth solution, as well as those 


where we obtain solutions with discontinuities. 


For simplicity consider a flat, inclined surface (z = x). 
The characteristics will be straight Lines in this plane, 
parallel to the x-z plane. Modifying the reflectivity of the 
surface to be increasingly darker with increasing x, we 
obtain a new solution which contains characteristics, again 


parallel to the x-z plane, but curving toward large z for 
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large x. This solution is smooth and contains no indication 


of an error. 


If now we apply Instead a surface coating which is normal for 
positive y and darker for negative y, we obtain a solutton in 
which the inconsfstency is apparent, The characteristics in 
the solution for negative y are more Inclined than those’ for 


positive y, anda discontinuity exists at y = 0. 


Using this kind of approach one could determine which kind of 
surface markings are noticeable by an observer (i.e. lead to 
inconsistencies in the solutton) and those which merely alter 


the apparent shape. 


2.15 WHAT ARE LIKELY SOURCE DISTRIBUTIONS? 


Since the complexity of the algorithm presented here 
increases with the complex ity of the Light-source 
distribution and stnce we only know how to bridge shadows 
cast by one source, it fs important to know which Light- 
source distributfons occur [n practice. First one notes that 
the situations found difficult by humans are almost certainly 
going to give difficulties to our algorithm. For example, 


when two sources cast shadows (such as on a road lighted by 
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widely spaced streat-lamps) the shape of unfamiliar objects 
becomes difficult to ascertain because of the crossed 
shadows. If the incident intensity varies greatly from one 
image area to another (such as in a lightly wooded forest) 
the tangle of Lighted and dark areas makes perception more 
difficult. On the other hand one would expect “natural” 
conditions to be particularly easy. That is, one point 
source somewhat above the observer (the sun) combined with a 
very diffuse (almost uniform) source (the sky). The diffuse 
source will not throw sharp shadows of its own, The absence 
of either of the two sources makes vision only slightly more 


qdifficult. 


2.75.7 RELEVANCE TO PHOTOGRAPHY AND GRAPHICS: 


One would expect photographers to have something to 
contribute to this subject and jfntroductory booklets on 
artificial light photography confirm the above conclusions, 
The beginner is advised to use anunber of lights with 
different characteristics as follows (Phrases of inexact 


meaning will be placed in quotes): 


1. The main light - The ideal main tight ts a targe spot 


light approximating the effect of the sun. It is usually 
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placed 45 degrees above and 45 degrees to the side of the 
subject. Its purpose is to establish the “form of the 
subject” and fix the ratio of lighted to dark areas. The 
exact ratio is not important but the position of the 
source should result in good shading (which increases as 
the source is moved further from the camera) without too 
much shadow area (in which detail ts more difficult to 


perceive). 


2. The fill-in light (Cor axial light) - Its purpose fs to 
lighten slightly the shadows cast by the main Light and 
approximates the effect of the sky. It fs placed near the 
camera to prevent it from casting its own shadows and to 
simulate the effect of uniform lighting (see an earlier 
discussion of uniform illuminatton, section 2.4.6). The 
appearance of shadows within shadows is considered 
extremely “ugly” and should be avoided since it makes the 
picture more difficult to interpret. The ratio of fill-In 
light Intensity to main light intensity ts usually chosen 


to be about | to 3, 


In addition a number of small sources may be used for extra 


effects: 
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3, The accent tight - Its purpose is to enliven the rendering 


by adding highlights and “sparkle’, It should be a small 


collimated source which can be directed to ftluminate 


small sections of the subject. It is placed behind and to 


the stde of the subject so that it cannot cast shadows of 


its own. This tight can add catchlights (specular 


reflections such as on eyes. or metal objects) and bright 


outlines (particularly on hair), 


4. The background light - Its purpose fis to “separate” the 


subject from the background, It illuminates the 


background only, such that the intensity reflected by the 


subject will nowhere match that of the background. This 


ensures that 


edge between 


Other hints are 
the main-light 
crossing edges 

light parts of 
It is 


formulated in 


interesting to 


the two can be easily “separated” - i.e. the 


them will be vistbte. 


that too many lights spoil the effect, having 


at the camera creates a “flat” image, shadows 


on the subject are to be avofded and that 


the image draw the attention of the viewer. 


note how much of what is vaguely 


these introductions to photography can be 


understood from the pofnt of view of shading, 
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2.16 DETERMINING SHAPE FROM TEXTURE GRADIENTS: 


A problem related to that of determining shape using shading 
is that of determining shape from the depth-cue of texture 
gradients. A textured surface will produce an Image’ in 
which the texture is distorted in a way reflecting both the 
direction and and the amount of the inclination of the 
surface, An image of a tilted surface with a random dot- 
pattern for example will be compressed in one direction (the 
average distance between dots is decreased) by an amount 
proportional to the inclination of the surface. Both 
direction and magnitude of the gradient can thus be 


determined - except for a two-way ambigulty. 


In practice it may not always be easy to determine = such 
texture gradients reliably because of tow resolution of the 
imaging device and scatter, causing a reduction In contrast. 
Some simple textures may be handled by simple counting or 
distance measurements as suggested above, while more 
complicated textures (e.g. a plastered wall) will need more 
sophisticated techniques, such as two-dimenstfonal correlation 
(best obtained using the fast-fourier-transform), Some 
experimentation with this technique showed promise, but did 


not supply very reliable gradients and the method was slow. 
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The next problem is how to obtain the shape from the texture 
gredients. Startins at sore point (whose distance from the 
camera we assure known), we use sore external knowledze to 
resolve the tworway ambiguity. Ve can now take a small step 
in any direction and find the gradient at this new point. 
Continuing in this way we trace out some curve on the surface 
of the object (somewhat analasous to the characteristics in 
the shape-from-shading method, except that here the curve’ fs 
quite arbitrary). 


Let s be the arc-distance along the curve, z#, the distance to 


° 
the initial point, and p and aq the components” of the 
gradient, then: 

S 


z(s) = z + Cigd)” a cs 


If one takes small enough steps, one can continue to resolve 
the ambiguity at each step by using the assumption of 
smoothness. This can be done until we meet a point where 
the gradient is zero. To continue past such a point would 
requlre some external knowledge to again resolve the two-way 
ambiguity. An aggregation of points with zero inclination 


can form an ambiguity edge which cannot be crossed. 
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Clearly we can reach a given point through many paths’ from 
the initial point. This allows us some error checking, but 
there certainly are better ways of making use of the excess 
information. For that is what we have, since we know from 
the solution to the shape-fromshading that only one value is 
required at each point for the determination of the shape, 
while we here have two (the components of the gradient). 
Most commonty when faced with such an excess of information 
on can make use of some least-squares technique to fmprove 
the accuracy, Perhaps a relaxation method on a grid would be 


useful (The grid need not be rectangular), 
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3, PRACTICAL APPLICATION: 


3.7] THE SCANNING ELECTRON MICROSCOPE: 


This chapter deals with a few practical applications in which 


the equattons simplify considerably, 


3.1.1 DESCRIPTION OF THE SCANNING ELECTRON MICROSCOPE: 


This device uses an electron beam which is focused and 
deflected much like the beam of a cathode ray tube and 
impinges on a specimen in an evacuated chamber [11], The 
narrow ray penetrates into the specimen for some distance, 
creating secondary electrons along its path (a small number 
of electrons are reflected at the surface). The depth of 
penetration, the spread and the number of secondary electrons 
are all functions of the materfal of that portion of the 
specimen. The number of secondary electrons which reach the 
vacuum through the surface will depend strongly on the 
inclination of the surface w.r.t. the beam, being least when 


it ts perpendicular, 


These relatively slow secondary electrons are then attracted 


by a positively charged grid and impinge on a phospor-coated 


ELECTRON 
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Figure 21: Sketch of a scanning electron microscope. 
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Figure 22: Detail of electron beam impinging on specimen. 
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photomultiplier. In this way a current is generated 
proportional to the number of secondary electrons escaping 
the specimen. There are other modes of operation which do 
not however interest us here. The output {fs used to modulate 
the intensity of the beam in a cathode ray tube while both 
beams are scanned synchronously In a T.V. like raster. The 
image created exhibits shading and is remarkably easy to 
Interpret topgraphically. This is quite unlike the normal 
use of optical or transmission electron microscopes which 


portray density and thickness, 


The magnification is eastly Increased by decreasing the 
deflection in the microscope, The resolution is poor 
compared to the transmission electron microscope because of 
the spread of the beam as It enters the specimen, but the 
depth of field is much better than that of an_ optical 
microscope because of the very narrow beam (extremely high f- 
number). The higher field gradient on edges causes these to 
be outlined more brightly. This artifact, while appealing to 
people, may be a problem in the implementation of a computer 


algorithm for finding the shape. 


Often the final analysis does not involve exact determination 
of the shape or two stereo-Images can be used, but there 


propably are also important cases where the shape must be 
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determined and the stereoscopic method fs not applicable. 
This may be because at the magnification used the specimen 
appears smooth without significant surface detail or because 
it is difficult to line up the second image. Since the 
equations for this case turn out to be so simple it should be 
rewarding to tie a scanning electron microscope directly into 


a small computer. 


3.1.2 EQUATIONS FOR THE SCANNING ELECTRON MICROSCOPE: 


A tittle thought shows that this is analogous to the case 
where the source is at the camera (or equivalently we have 
uniform illumination); for one thing, no shadows appear. 
Next we note that at all but the lowest magnifications the 
projection fs near-orthogonal. Because of these two effects 


the five 0.D.E.%s simplify considerably: 


x= Fe, y=, 2 pFp + GFa 
p = ~FL - pF, and q = ary - aF, 
Now F, = A o,I, and F, = -b; 
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l= nein S 1 /n (where n =(-p,~q,1)) 
In= (/n)(2 - 1A) = (B/n) - (1 /n 30K 


I, = (I/n?)p and Iq = (1/n*)q 


Hence: k = Fp = (A ¢,/n3)p, ¥ = Fy = (A dg, /n)q 
(A $,/n?)(p* * q*) 


-b 


Me 
i) 


x and qe a ore 


Te 
ii] 


If ds # 0 everywhere, we can change to a new measure s along 
the characteristic by multiply ing ail equations by 


dN = nPs(A x) and we get: 


K= Pp, Y=a, =p +q° 
p= by (n° /CA b.)), a = b(n /CA O,)) 


This extremely simple case thus has characteristics which are 
curves of steepest descent (or ascent). Also note that the 
equation for 2z does not couple back into the system of 
equations (due to the orthogonal projection) thus increasing 
accuracy. The equations happen to be very similar to the 
efkonel equations for the paths of light-rays in refractive 
media. It may be possible to find ready-made solutions to 


some special cases by using this analogy. 
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We assumed that oy #0; this is equivalent to assuming that 
an inverse exists which atltows us to find I from a 


measurement of the image intensity: 


I 
(1- Cx) /(2WOx)) 
(1/2)*(p~ + q™) 


W(o(1,1,1)) 
Let 30x) 
Then Z(O(1,1,1)) 


So we can find at each point the magnitude, but not the 
direction of the local gradient. This is very different from 
the method of determining shape from texture gradients 
(section 2.16), where we can locally determine the gradient 


except for a two-way ambiguity. 


3.1.3 AMBIGUITIES AND AMBIGUITY EDGES: 


This is an easy enough example to. study amb Igufties, 


Consider the two surfaces: 
3 
z5=z+¢ x , 222+ [x] 
Clearly they cannot be distinguished from monocular views 


since their gradient magnitudes are tdentical: t.e. they 


produce identical Intensity distributions in the image. This 
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Figure 23: A locally determined ambiguity edge. 
f=1/ (x +y -1)" 


ie 


an 


Figure 24: A globally determined ambiguity edge. 
f= 1/(14x? + (y-1)?)  1/(14x? + (y+1)*) 
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manifests itself in a slowing down of the characteristics as 
they approach the Line x = 0 (alternatively SS oo ). They 
cannot cross this line aggregation of singular points, Note 
that the characteristics approach this tine at right angles 
and that the edge ts determined locally, each point on ft 


being a singular point. 


A second kind of ambiguity edge can occur parallel to 
characteristics, separating those which can be reached from 
one singular point from those reachable only from another. 
This kind of edge is not tocally determined, since a change 
in the surface is possible which removes one of the singular 
points and makes all the characteristics accessible from the 
other. This can be done without altering an area near two 


given points previously separated by an ambiguity edge. 


Both types of ambiguity edges occur in the general case but 
are not so easily studied there. They divide the image into 
regions within each of which a_ solutton can be obtained, 
Typically most such regions will have one singular point from 
which one may obtain initial conditions (provided one makes a 
decision about whether the surface jis concave or convex and 


knows the distance to the singular point), 
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3.2 LUNAR TOPOGRAPHY: 


3.2.1] INTRODUCTION TO LUNAR TOPOGRAPHY: 


The other very interesting simplificatton to the general 
shape from shading equations occurs when we introduce’ the 
special reflectivity function which applies to the material 
in the maria of the moon. This in fact was the first shape 
from shading problem solved both theoretically and in an 
operating algorithm [4]. Using the special reflectivity 
function and the fact that the sun is a distant source, it is 
possible (but very tedious) to show that the equations 
simplify so that the base characteristics (i.e. the 
projection of the characteristics on the image plane) become 
straight tines radiating from the zero-phase point. This 
point corresponds to g = 0 and is directly opposite the = sun 
as seen from the camera. Actually this is true only when the 
sun is located at negative 2, for positive z (that Is in 
front of the camera), the relevant point ts the phase 


point, directly fn the sun. 
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function. In most of these formulaes we find that for a 


given G, b is constant for constant I/E. The tines of 


constant I/E are meridians, 


At full moon, when G = I we find that the whole face has 
constant Luminosity. This is quite unlike the effect on a 
sphere coated witha typical matt paint where the image 


intensity would vary as: 


/1 - (r/R)” 


Where R is the radius of the image andr the distance’ from 
the centre of the image. The full moon thus has the same 
appearance as a flat disc if one is used to objects with 
normal matt surfaces, This may explain the flat appearance 


of the full moon. 


3.2.3 DERIVATION OF THE SOLUTION FOR LUNAR TOPOGRAPHY : 
3.2.3.1 THE BASE CHARACTERISTICS: 


In the case of pictures taken of the lunar surface from 


nearby (e.g. from orbit) we have the following: 
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1. Distant source (the moon subtends an angle of about .03 


milli-radians at the sun). 


2. Near point source (the sun subtends an angle of about 10 


mitti-radians at the moon), 


3. Camera at the origin. 


4. The reflectivity function is constant for constant I/E. 


This is a property of the material of the maria of the 


moon which has heen known for some time, 


we have (using results obtained in subsection 2,3.4): 


a 


I, = 0 In = (1/n)(r,- 1) 
Ep = (I/r)(h- — F) EE, = (s/n) (> EF) 
Gr= (/r)(P- Gr) Gye 0 


Where Fe is a unit vector in the direction from the sun to 


the moon. 


If I and E depend on some parameter s, while I/E is constant: 


Since (1,E,G) 
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is constant for constant I/E: 


- Le b. E, = 0 and therefore: 


Le, +E be = 0 


If I and E depend on soma parameter k: 


dily + 6, Ey 


b Cle = (IER 6. (EET = IE) 


bs - o,£ caus 9 cag! Pe 


ee Oa 


be Fe = Ebay (Ely ~ TE) 


Using some of our previous results we find: 


+ 
nt 


-(I/r)(n - E 4) 


n = (E/n)F, = (EI/n)R = CI/n)F + CEL/n)a 


= (1/n) (CRAP - CF ADT] = Cl/n) (hx Fx 


And since A, = 0: 


We will ignore 


tT 
tt 
> 
~~ 
t 
o 
ly 


n 
Us 
i" 
so 
> 


FL for now, mainly because it has an ugly 


Looking expanston. 
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Ix, Cz-ay )-x(z,~Qy,) ,¥, Cz-px )-y Cz Q-PxXo ), 20x pty 09g )-z, (xptyq) ] 


Let X = -AE o,, 


“aw 
where (x,,¥, +20) = Foe 
n*rr, 


Note that A is a constant in this case, 


Xo y x Yo 

Fp = X}— (1 - a-) - -G = a4 
Zo Zz Zz z 
Yo x y Xo 

Fa = X (1 - p=) - -(1 - p) 
Z 5 z z z 


Now Looking back at the five 0.D.E.°s: 


+ ay 


xe 
" 
n 
7 
~ 
We 
tl 
mn 
- 
. 
™ 
ft 
ae) 
mT 
7 
+ 
co) 
n 
~~ 
i 
v 
x 


oe 
t 
i 
x7 
ae) 
” 
y 
. 
Qe 
iL 
i] 
nn 
~< 
t} 
a 
we 


Again we can decide to tgnore p and q for the time being, and 
attempt to determine the behavior of the characterlstics, 
Our aim is to show that their projections fn the image plane 
are straight lines independent of the scene. The behavior of 
y against x Is of Little help and we next took at the 


projections in the Image plane: 
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x’ = x (f/z) and y’ = y (f/z) 
O/#) = (1/22) haexd) = (1/2) hanx (pxtay)] 


x e X oe 
(1 - -p)x - ay 
Zz z — 


y X Xx y 
p= gs) =O eps = 9-2 
z z z z 


x” Xx, x x y 
- = -{¢ - -)(1 - p- - q-) 
f z\ z, z zZ F 


yo Xfy, y x oy 
~ = -|(—- -)(1 - p- - aq-) 
f z\ Zo F Zz z 


Now if the surface is not tangent to the ray from the camera: 


r 
mot 
m | 
o L-¥ 
v~_ 

' 
mo 1X 


Stmitarly: 


E # 0 i.e. ren 7 0 and therefore: 


x y 
(1 - p- =~ q-) # 0 
Zz z 


If in addition bax #0, A # 0 and z # 0, then we can divide 


the two equations: 


° Se 
° 

y Zo z 
ee = 

x Xo x 
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This first-order ordinary differential equation for the base 


characteristics is separable: 


-dy’ -dx* 
Ver Ve Xe Xx 
Zo f Zo f 
Solving this we obtain: 
Y =Y ea 
log (—— - —) = log(— - —) + log (c) 
Zo f Zo f 


Let the arbitrary constant c be tan(t): 


1 ee i ee ee 
sin(t) [z,  * |” cost) lz, f 


Thus the projections of the characteristics are straight 


Lines in the image plane emanating from the point: 


Bs. see 
Z| ze 


If the sun ts behind the plane of the image (z,> 0 - as 
would usually be the case for reasonable illumination and 
avoidance of extraneous light entering the lens) this point 
is called the zero-phase point, since {it corresponds to the 
point In the scene which is directly opposite the sun as seen 
from the camera and hence g = 0, Because of the special 
properties of the reflectivity function of the maria of the 


moon intensity variations in this region are entirely due to 
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non-uniform surface properties rather than shape. It is for 
this reason that this point ts not usuatly tncluded In the 
image but lies somewhat outside it In the image-plane, This 
will prove unfortunate later on when we have to invent 


initial conditions. 


If the sun is in front of the image plane (z, <¢ 0), the 
special point ts the TT phase point, where the image of the 


sun would appear in the image-plane. 


So the obvious simplification to the equations which would 
arise if we let xo. = y, = 0 cannot be exploited since we do 


not wish to ortfent the camera in this way. 


We would like to arrange for s, the parameter that varies 
along each characteristic, to correspond to arc-length. This 
can be achieved by multiplying each of the five 0.D.F.%s by 


» , where: 


ane Fe 1 
A= X's x y 
Closepe nga) 
z z 
x fk x) Lo xo) 1 
f Zz Zz] S z f } s 


Then by choosing constants suitably: 
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x x X, 
- = = + s cos(t) 
Ff Ze 
@é 
y y y 


= —+s5 sin(t) 
z f Zo 


Thus s gives arc length along the characteristics while the 
value of t selects a particular characteristic. 
3.2.3.2 THE INTEGRAL FOR 2z: 


We next turn to £ which we would like to find without solving 


e e 
the messy equations for p and q. 


e e ® Xo Xx Ye y 
Z = px + py = XX Ip(— - -) + q(—-o- -) 
Zo Zz Zo Zz 


= XX s (p cos(t) + q sin(t) ) 


This is a good place to Introduce some abbreviations of 


commonly occuring dot-products: 


Xe Ye 
(—,—) . (cos(t), sin (t) ) 
Zq Zo 


r- 
! 


M= (p,q). (—,—) 
Zo Zo 
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N = (cos(t), sin(t) ) . (p, @) 


Note that L is predetermined (i.e. independent of the image) 


and that L and M tend to 0 if the camera is pointed directly 


away from the sun (i.e. Xo = Yo = 0) 
xX y 
(1 - p- - q-) = (1 -M-sN) 
z z 
Zz -] 
X = -* 


s (1-M-sN) 


‘ ~zN 
and so: Zz = 
(1-M=sN ) 
We now attempt to express this in. terms of measureable and 
calculable quantities (s.a. G, I/E, s and t). Since br% # 0 
and differentiable it must be monotonic and hence have an 
inverse. That Is, given b/A we will be able to calculate I/E 


(G is known at each point). 


te = (%o-¥,e%) and rr, = [xr + yr+z2 


Xe 


1 ie (x,y,z) = 2( 


+ s cos(t), 
Zo Zo 


r 
Let Qe / s* + 2s + (2° 
- Zo 


+ s sin(t), 1) 
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VY (b/A,G) = I/E 


b r, | 
(—,G)#(—)- - 1 
" A Zz, Q 


The usual tables for b in the case of the marta of the moon 


ws 
" 
nw tm 


however are not usually given in terms of I/E and G, but 

rather & and eg. Where : 

G-I/E 

tan(d&) =——— 
1-G* 

« is the projection of the emittance angle on the phase- 


angle plane. 


z T z, (I-M) 
G- I/E = —x—- Q —————_. 
r, Q rm (1-M-sN) 

z, / r. CHEM? 6 
= oe Pea a) 6k.. Se 
r, Q Zo (1-M-sN) 

z § N es 
2S St 1 (Stl) eee 
lo Q (1-M-sN) 
ae 
1-G = 7 - 
ro Q 


eee 
De f ine P = sgn(z,) /(—) ob 
Zo 
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G-I/E 1 N - 
——— x = Sa a a Sey 
Fi=G* P \(1-M-sN) 


3.2.3.3 THE INTEGRAL FOR r: 


So far we have been working in the coordinates x*, y° and 2. 


The final result looks neater if we use r Instead of 2. 


r = z2Q 


z2Q + z(s+L)/Q 


me 
Li} 


zN (s +L) 
co Q + 2 
(1 -M-sN) Q 


r N 2 
el ay (s+L) 
Q \(l-M-sN) 


~(rP/Q™) tan() 


s 
u 


Written out more fully, we have: 
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PWG cs, ae 
sen(z,) /(—)° -L* 
Zo 


= _ tan (XK) 
r. . ‘ 

s +2sl+(-——) 
Zo 


a | 7. 


The numerator is a fixed quantity for each characteristic, 
the denominator varies along each characteristic (but is 
independent of the scene), while tan(&%) is obtained from the 
measurement of b/A and the known G (using the function W). 
The given ordinary differential equation for r has the simple 


solution: 


s 
tan (*) 
-Pp ds 
r(s) ° Q” 
=e 
r(0) 
Xo Yo 
where L = — cos(t) + — sin(t) 
Zo Zo 


y Fy ro 2 = 
and Q=/s +2sl+ (—) , P = sgn(z,)/(—) - L 
Zo Z 


r(0) is the distance to the point from where the integration 


was started. 


To sum up: as one advances’ along each characteristic in 


turn, one calculates G, measures b/A and uses WV to obtain 
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tan(%), which is then used in the evaluation of the above 
integral. The process is much simpler than the general shape 
from shading algorithm in that the base characteristics are 
predetermined straight lines and only an integral needs to be 
evaluated, It is possible to write the above’ result ina 
slightly more elegant form, which is the one derived by T. 


Rindfleisch (for z,) > 0): 


' $ 
cate eRe eT aae (r*. 2) tan(%) ds* 
r(P) fléxi,f 4B 7 
= e o. = 
r(Py? 
Whe re s° = fs 

N,= mee” = -(f/z2) (rere) 

f y sy x x xy x oY 
= -- zz,( (- - —), (2 - -), Hi - >) ) 

z Zz zo Zp = ZZq «ss: Zot 


y x, 
-~fz,s(-sin(t), cos(t), —cos(t) -—sin(t)) 


2. Zo 
2 x L 
- 0 “0 
N ow L « — cos” (t)+2 cos(t) sin(t)+— sin (t) 
Zo Z5Zo z3 
Ae aN, xe 
No = fz s/ i+—cos (t)-2 cos(t) sin(t)+ sin (t) 
Zz; Zz, z> 
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2xN, = -fz js (cos(t), sin(t), 0) 
[2x] = sflz,| 
1 No 
= a = Pp 
| 2xN | | ZN, | 
(P7380 = CRB & CreB) fr™ = Ce/r)™= 1/0” 


The two ways of writing the integral are thus equivalent. 


3.2.4 SOME COMMENTS ON THE INTEGRAL SOLUTION; 


The base characteristics are predetermined straight Lines 
(independent of the image). This makes for high accuracy 


and ease in planning a picture taking mission. 


Onty asingle integral needs to be evaluated, not five 


differential equations. 


The primary input is the intensity, not its gradients, 


again making for high accuracy. 


Although, as usual, the reflected light-intensity does not 
give a unique normal, it does determine the slope 


component in the direction of the characteristic. 


Page 110 


J. van Diggelten [2] first noted a special case of this 
when he solved the lunar topography problem for the 
special case of an area near the terminator (Line 
separating sunlit from dark areas). The characteristics 
are such that the slope along them can be determined 
locally. The slope at right angles to the 


characteristics cannot be determined locally. 


5. Although T. Rindfleisch [4] did not mentton it in his 
paper it is very easy to bridge shadows since each light- 
ray Lies in a sun-camera~characteristic plane. Its image 
can thus be traced on the base characteristic until we 
again meet a lighted area. One need not even make special 
provisions for this, but just use tan(%) for grazing 


incidence (intensity = 0) in the shaded section. 


3.3 APPLICATION TO OBJECTS BOUNDED BY PLANE SURFACES: 


Since a great deal of image processing these days is applied 
to images of polyhedra one might enquire how one could apply 
this method to such objects. First we note that the main 
features of these objects, the joints (angular edges on an 
object) and edges (where one object occludes another), are a 


stumbling block to the application of our method developed so 
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far. Since we already know that the areas of more or less 
constant reflectivity are plane faces there is Little point 
in exploring them, So acompletely different approach is 
indicated. Firstly we might check whether a parsing of the 
scene obtained by some other method is correct In the sense 
that the intensities observed for the faces correspond to 
their inclinations (the information on the intensity of the 
faces is usally discarded), It is not clear however what one 
might do if this test fails. Furthermore, the programs which 
reduce the image to a ltne-drawing and then decide whieh 
faces belong to each object cannot really determine the 
inclinations of the varlous faces without additional 


assumptions (orthogonality for example), 


One can however find the normals to each face directly from 
the known slopes of the projection of the joints In the image 
and the measured reflectivities. One must know which Lines 
In the image are true joints (between two faces belonging to 
the same polyhedron) and which are fortuitous (edges between 
faces of different objects). For each normal we need = two 
values, Each intensity gives us one non-linear equation and 
each slope of an image of a joint gives us another, The 


equation from the intensity is of course: 
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AC(r) bere,@) - b(r*) = 0 


Where we know that I, E and G are functions of p and q. There 
wilt be one such equation for each face. 


Where two faces with normals n, and n, intersect, they forma 


\ 
joint which will be seen in the image. Suppose two points on 
this image are A and B. Then a vector perpendicular to the 


plane through the joint and the camera is Axes, We also know 


that the joint must be parallel tonmxn,. But AxB must be 


~ wf 


perpendicular to the joint hence: 


IMAGE 


Figure 25: Projection of a joint on a polyhedron on 
the image. 
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Each joint contributes one such equation. Next we determine 
how many faces must intersect before we have enough 
information for a solution. Two faces intersecting give us 


two equations from the intensities and one from the stope of 
the image of the joint, while we need four unknowns, An 
infinity of solutions thus exist, With three faces a 
solution is possible since we have six equations in_ six 
unknowns. Because of the non-linearity of the eauations 
more than one solution might exist. With a larger number of 
faces we always have at least enough information fora 
solution and at times have more equations than unknowns which 
may remove some of the remaining ambiguities and improve’ the 
accuracy. In this way too it may be possible to discover 


which joints are really between faces of the same object. 


3.4 FACIAL MAKE-UP: 


When a surface whose photometric properties are taken to be 
uniform is treated so as to change these properties in some 
areas, the apparent shape is changed. This of course is one 
of the uses of make-up. The shape of a face for example can 
be made to conform more closely to what a person thinks’ is 
currently considered “Ideal”. This fs achieved by making 


some areas darker (causing them to appear steeper) and others 


=. 
~ 
bed 
Ye 
by 
G 
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Illustration of the effect of facial make-up. 


Figure 26: 


-up. 


Illustration of the effect of facial make 


Figure 27: 
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lighter. Areas lightened usually include singular points and 
cause a change in the apparent skin darkness (a normalization 
effect) and will change the apparant shape in areas other 


than the singular points, 


These modifications can change the shape perceived when 
viewed under the right lighting conditions. The effect will 
change somewhat with orientation and may at times disappear 
when no reasonable shape would give rise to the shading 
observed, Because of a number of surface oils the skin has 
a specular component in its reflectivity, it fis also fairly 
translucent. Both of these effects are sometimes controlled 
with talcum powder, The removal of the specular components 


makes the surface appear more smooth. 
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4, EXPERIMENTAL RESULTS: 


4,1 A PROGRAM SOLVING THE CHARACTERISTICS SEQUENTIALLY: 


When the solution to the shape from shading problem had been 
found using the inconvenient coordinate system (x*, y*, z), a 
program was written which would solve the five O.D.FE.°s along 
one characteristic at a time (the equations used are not 
reproduced here). The input data was obtained from the 
image-dissector camera attached (at that time) to the PDP-6 
computer in the Artificial Intelligence Laboratory. This 
camera is a random access device: when given an x anda y 
coordinate it returns a number proportional to the intensity 
at that point in the image. The program first searches for a 
maximum in intensity, constructs a small spherical cap around 
it (to obtain an initial curve) and uses a standard numerical 
method (see subsection 4,171.2) to solve the set of five 


ordinary differential equations, 


The prime data required in this soltutton is the intensity 
gradient (x* and y” derivatives of the intensity), which is 
obtained from the intensities measured for a small raster of 
points near the current x* and y”. A linear function In x’ 


and y is fitted to this set of intensities; the 


coefficients of x” and y” are the desired gradients. The 
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size of the raster is chosen to correspond to the step-size 
used in the nutrerical solution method, so that successive 
rasters almost, but not quite, touch. In this way fair 
accuracy in the determination of the gradients is obtained 


without loss in resolution. 


If the teast-squares fit is bad, indicating that surface 
detail is being missed with the stepsize used, or that the 
characteristic is traversing an edge or joint, the solution 
for this characteristic is halted and the solution started 
for the next characteristic, Other reasons for terminating 
the solution are that the characteristic has left the field 
of view of the image-dissector or reached a very dark region, 
most likely a shadow or the background. When either the 
(catculated) incident or emittance angles become very small 
(indicating approach to an edge or shadow edge) or »~ very 
large (indicating approach to another singular point or an 


ambiguity edge) the solution will also be stopped. 


The data structure here is very simple; just a record of 
varlous values ( x*, y’, z, intensity, p*® anda*) for each 
point on each characteristic. The shape so determined can be 
displayed in perspective and stereo ona DEC 340 display. 
The characteristics appear as dashed lines - each dash 


representing a step in the integration (We chose the 
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parameter s so that each dash represents the same arc 
length). The output can be photographed from the display 


and plotted on a Calcomp plotter. 


4.1.1 AUXILIARY ROUTINES: 


A number of auxiliary routines needed to be written for this 
program. First the Incompatable Time Sharing System (ITS) 
on the PDP-6 does not support a FORTRAN style arithmetic 
language and all programming was done in assembly language 
(MIDAS). The large amount of arithmetic Involved, 
particularly with the inconvenient notation and coordinate 
system used at first, made it imperative to incorporate into 


the assembler the ability to handle arithmetic statements. 


Next we constructed a package of useful routines which 
handles floating point I/0, dynamic array allocation and easy 
generation of display tists for the DEC 340. It also 
includes routines for the standard arithmetic functions 
(SQRT, SIN, LOG etc.) and manipulation of vectors and 
matrices (multiplicat ton, addition, inverston etc.). 
Interrupts, user defined operations and command 
interpretation are dealt with as well, Some of the remaining 


routines will be briefly described in the next few sections, 
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4.1.1.1 STEREO PROJECTION AND OBJECT ROTATION: 


Since it ts important (particularly during the debugging 
phase) to be able to visualize the shape being calculated, 


stereoscopic output on the display fs provided, 


Let de be the separation of the eyes, f their distance from 
the display and dy the distance fromthe eyes to the origin 
of the coordinate system (usually chosen to be at the 
singular point from which the solutton was started), The 
coordinates of the left-eye and right-eye images of the point 


(x,y,2) are then (xoey) and (x/,y*) where: 


Ko e 7 (xtds/2)(F/(2tde)) Fds/2 
| 
y’ = y (f/(Cz+d,)) 


A pair of lenses is employed to focus on the surface of the 
display while converging on the apparent point (x,y,z). 
Obv lously one needs to know the scaling of the display in 


terms of dsiplay units per mn. 


One would like to be able to view the objects from various 
sides and perhaps even have some rotational motion to gain a 
greater perception of depth. To this end the object can be 


rotated around its origin (also offset and expanded tn size). 
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DISPLAY 
SURFACE 


Figure 28: Stereo-projection of an object point. 


ROLL Zz 


PITGH 


YAW 


Figure 29: Definition of Pitch, Yaw and Roll. 
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This was preferred over the more common method of allowing 


the eyes to be moved around in the object space. 


To obtain any orientation with one matrix multiplication, the 
three angles P (pitch), Y (yaw) and R (roll) are defined as 
rotations about the x, y and z axes respectively. They are 
applied in that order (order is important because rotations 


are not commutative). 


cos(R) -sin(R) 0 cos(Y) 0 sin CY) 1 0 0 
A =|sin(R) cos (R) 0 0 i 0 0 cos(P) -sin(P) 
0 0 1 -~sin(Y) 0 cos(Y) 0 sin(P) cos(P) 


Using the abbreviatton c for cosine ands for sine we have: 


ck cy cR sY sP - sP cP CR SY cP «+ sRi-sP 
A =| sR cY sR s¥Y sP + cR cP sR sY¥ cP - cR sP 
Sy cY sP cy-eP 
The various parameters controlling the object rotation and 
the projective mapping can either be preset or continously 
read in from anumber of potentiometers (connected to a 
multiplexor and an A/D convertor) controlled by the viewer. 
While one display list appears, the other fs being calculated 
using the latest set of parameters andwill in turn be 
displayed when completed. The parameters are also displayed, 


they are: 
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of the Llight-source is measured, as well as the distance from 
the front of the sphere to the entrance pupil, The imase 


dissector ts focused on the edse of the sohere. 


With the sphere temporarily illuminated from several sources, 
the program finds its exact position and size, as well as the 
difference in horizontal and vertical deflection sensitivity 
of the imape-dissector. It is now In a position to calculate 
the pcints tn the image which correspond to given incident 
and emittance angles. For aonumber of choices of both of 
these angles it then reads the intensity at a small raster of 
points (to reduce netse and the effect of pineholes in the 
photo-cathocde) near these positions and averages then, Since 
there are usually two places in the image with the same 
incident and emittance angle, a check on the data is 
available. The resultant table of values (usually 
normalized w.r.t. the brightest intensity) can be printed and 
the whole process repeated after moving the lLight-source to a 
new position for a new phase angle. The program accounts for 
such things as change in incident light intensity as the 


light-source gets moved around, 
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4.1.1.3 FINDING THE CALIBRATION SPHERE: 


This subsection and the next deal with details, needed In the 
program for measuring the reflectivity function, which may 


not be of general interest. 


For good accuracy we first need to know the exit pupil to 
image plane distance (the focal length is given). This would 
be easy If one could focus on the front of the sphere. It 
turns out that a simple approximation will work ina few 
iterations. At each step one recalculates the estimated 
distance to the edge of the sphere, the estimated exit pupil 
to image plane distance and the estimated radius of the 
sphere, using the previous estimates and the measured radius 


of the image. 

Next we need to find the exact center and radius of the 
sphere from its image coordinates and the known distance to 
its front. First constder hortfizontal coordinates only. 


x, = f tan(a), x, = f tan(atb) and xz = f tan(at2b) 


We are given x, and x3 and wish to calculate x,, which can be 


done after expanding the tangents. 
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Figure 30: Determining the exact position of the calibration 
sphere. 


Figure 31: Finding points for given incident and emittance 
angles. 
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tan(2b) = (xg-x,) /F=x 3X1) 


tan(b) = Si + tan*(?b) - T)/tan (2b) 


x,/f = (x, + f tan(b))/(f - x ,tan(b)) 


The same formulae are then used to find the vertical position 


y,. Finally we need to find z,: 


r= Retan(b) (/1 + tan*(h) + tan(b)) 


= Rg f/f ee 


zB (UR, #6) /Rs 


® 


Nm 
iT} 


4.1.1.4 FINDING POINTS FOR GIVEN i AND e: 


Clearly the points for given incident angle lie on a circle 


(on the surface of the sphere), Similarly for points with a 
given emittance angle, These two circles may intersect in 
two, one or no points, One can find this intersection by 


first finding the Line along which the planes containing 
these circles intersect, Applying the sine and cosine laws 


we get: 
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Let I = cos(i) as usual and D= ysl 
D/sinQ™i) = r/sin(a) b = isa 
r cos(b) = r(cos(i) cos (a) + sini) sin fa)) 
=(r/0) cr /o* - Pas 1. ee eS 9 
d = r cos(b) 
a Me fad Ve 


The equation of the plane In which the circle of points with 


given incident angle i lies is: 
nw a 
NON go ONS OE Say (where vi = (x,y,z) ) 


One can find a similar equation for the plane in which the 
circle of point with given emittance angle e lies. The 
introduction of an arbitrary third plane allows us to find 
one point von the intersection of the first two. The Line 
of intersection of the first two planes must be parallel to 
the cross~product of their normals (let them be v,, and vy, ). 


So the equation of the line we are looking for is: 


The points we are trying to find must also lie on the sphere, 


i.e.: 


Page 128 


2 2 
(v- v,) =r 
~ ~ 
a L 
CV Ry Ao Net rat 
LL 
k*v we ot 2 KV, nm ve) + (vo ev) = rl = 0 
{ { al nl pon € 


The above equation may have no solutton for k, In which case 
no point exists for the given incident and emittance angle. 
Otherwise we can use the two solutions and substitute back to 
obtain the desired coordinates which are then transformed 


into image coordinates, 


4.1.1.5 SOME REFLECTIVITY FUNCTIONS: 


The first paint investigated was a matt white paint 
consisting of particles of Si02 and TiO, suspended in a 
transparent base. Very roughly one finds that the 
reflectivity function behaves like cos(i) fora given g. 
After playing with polynomial fits for a while, the following 
fairly accurate formula was found by aoprecess of little 
interest here: 
(1+#6)(2+G) 1+21EG-(I™+E 4G") 
4(1,£,6) = [1 + 
6 16 (1-G) 
Note the appearance of the discriminant discussed in an 


earlier section (2.7.3). The symbolic manipulatfon program 
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Figure 32: 


Table of reflectivity (for a white matt paint) versus 
I = cos(i) and £& = cos(e) for G= cos(g) = 0.8! . The 
intervals chosen correspond to constant size steps [in the 
angles. Note the blank areas for combinations of angles 


which cannot forma spherical triangle (see section 2.1.3) 
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MATHLAB [13] was used to find the derivatives bo be and be 
needed for the shape-from-shading program. For “reasonable” 
angles the above formula is about 5% accurate, becoming worse 
for extreme angles. The repeatability of this measurement 
was disappointingly low, depending on the depth of the paint 
coat andthe conditions of its application. Much of the 
investigation of the behavior of the image-dissector was the 
result of efforts to trace the remaining causes of 


Inaccuracy, 


Some other paints and an egrshell showed a matt component 
similar to the above, plus a very strong specular component 
(which is small except near the point for which i = e and T— + 
e = g). This component ts very sensitive to small changes in 
the surface properties such as can be brought about by 


handling the object. 


The image of a convex object with such a surface will usually 
have two local maxima in intensity. One of these will be 
broad (corresponding to the matt component), the other narrow 
and bright (corresponding to the specular component). These 
may be distinguished by a computer program on the basis of 
just these properties. It would then be possible to start a 
solution from the matt maximum (which is not ae global 


maximum) rather than the specular maximum. This might be a 
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good idea because of the increased accuracy (for one thing 
the normatization of image intensities would be more 


accurate). 


For the nose~recognition program, a plaster nose was used 
Initially, coated with the matt paint described above. This 


of course was not suitable for the final experiments. The 


) 

4 
1. A matt paint. 
2. Lambertian reflec. 
3. Skin on nose. 

3 
ro) I 
2 4 


Figure 33: Comparison of some reflectivity functions. 


restricted lighting conditions described later were chosen 
partly to avoid having to find the full-flegded function o of 
three angles for. skin. Since no true sphere covered with 
skin is available, measurements were taken of the shape of a 
real nose and intensities in an image (of a transparency) 
used to est imate p(LL1), In this way the non-linearities 
of the sHebaetsph ic process (and they were great) did not 


have to be determined separately. The properties of skin are 
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of course not very uniform and also vary fron person to 
person, so no great effort toward accuracy was made, Skin 
has a highly variable specular component, sO any 
normalization had to be done not w.r.t. the brightest point, 
but one nearby, The resultant table of (1) versus I lies 
somewhat below the one obtained from the matt paint under the 


same lighting conditions, 


4.1.1.6 PROPERTIFS OF THE IMAGE-DISSECTOR: 


In an attempt to track down poor resutts fn the first try at 
finding reftectivity functions accurately, the image- 
dissector was investigated in some detail [9]. Amongst 
problems found were: 

1. Unequal deflection sensitivity in horizontal and 


vertical directtons (differed by 12%). 


2. Twist of image varying with distance fron center of 


field of view, 


3. Poor resolution (3 Line-pairs/mm - radius of tube 50 


mm). 
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Figure 34: Geometric distortion in image-dissector for a 
triangular raster of points covering the photo- 


cathode. (The arrows are exaggerated 3 times.) 
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4, Pinholes in the photo~cathode (about 20 of up to 0,5 


mm in size). 


Non-uniform sensitivity of the photo~cathode (varies 


Ww 
. 


more than 30%). 


6. Fairly long settling time of the deflection coils, 


7. A large amount of scatter, which reduces the 
contrast by almost one~half and causes intensittes 
measured on the imase of a uniform square on a dark 
background to vary by 20%, depending on how close to 


the edge the measurement is taken, 


Some of these difficulties are inherent in the state-of-the- 
art of these devices, others were repaired, In any case, I{t 
was possible now to think about how to Improve the program to 


be more insensitive to these shortcomings, 


The program for finding reflectivity functions using spheres 
as calibration objects was sensitive to the (at that time) 
severe deflection inaccuracies, since the emittance angle 
varies rapidly near the edge of the sphere (this effect could 
be reduced with a parabolic test-object), A calibration 


table was created by another program in which are’ recorded 
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the image-dissector coordinates of a rectangular raster of 
equally spaced points on the photo-cathode, Also recorded is 
the sensitivity of the photo~cathode at each point. A simple 
interpolation routine can then be applied to coordinates sent 
to the image-dissector to counteract the distortion and, 
similarly, the intensity values returned can be corrected, A 
more conventent triangular raster of pofnts covering the 
whole photo-cathode was later established, Adjustments to 
the image-dissector eventually reduced the distortion by a 
significant factor and use of this table was no longer vital, 


although it did improve accuracy, 


4.1.2 NUMERICAL METHODS FOR SOLVING THE 0.D,.E.°S: 


The five 0.D.E%s were at first solved using a wel! known 
Runge-Kutta method [7, page 212]. The idea Is that ata 
given point we can calculate the derivatives of the five 
variables (x%, y”, z, p* and a’) w.er.t, the parameter s, 
Using these we take a half-step forward (increment s by h/2 
and calculate new values for x%, y’, 2, p” and q° as though 
derivatives higher than the first where zero). We then 
calculate the derivatives at this new point and take the same 
step (now using the new derivatives, which witt differ 


slightly from the previous ones). We next take a full step 
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Figure 36: Same solution viewed after a slight rotation. 
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(increment s by h) The final full step is taken using a 
weighted average of the four derivatives found in this way. 


Written out fn symbols this becomes: 


Let h be the step~size (for the parameter s) 

And Y = (x*%,y%,z ,p°,q") 

Also let the equations for the derivatives be: 
x = F(s,Y) 

(In our case F is actually Independent of s) 


th 


Denote Y(s,) by Y, then at the n step: 


Ky o h FCs. Yu) 

Ki = h Fs +h/2, Y,#K,/2) 

Ka = h F(s,th/2, Yn+K,/2) 

Ky = h F(s,+h, ¥, *Ks) 

Yar = YQ + 1/6) (K, #2K 42K 54K y) 


This method is easy to start (requires no previous values of 
XY) and stable, but requires four time-consuming evaluation of 
the derivatives per step. For this reason vartous predictor- 
modifier-corrector methods [7, page 194] were tried and the 


simplest was found to give adequate accuracy: 
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Pra = Yn + 2h Fls,, Y¥,) 
Maa = Raa 7 €4/500P, - Cy) 

Cue = Yn + Ch/2)CFCs,, Mya) + FOsy, Y,)) 
Manat = Oya t+ 1/5) CP ua - yar) 


P, M_= and Cc are the predictor, modifier and corrector 
respectively, This method is stable and requires only two 
derivative evaluations per step, but Is not self-starting. 
The Runge-Kutta method was retained for the first step in the 
integratton. Stability and accuracy were not. serfous 


concerns since the notse in the data Input contrtbutes far 


more to errors fn the solution. 


4.1.3 ACCURACY OBTAINABLE: 


Under optimal conditions (using the methods described to 
cancel out most of the distortton and non-uniformity of 
photo-cathode sensitivity) the program was allowed to scan a 
sphere of 100 om radius, A sphere was then fitted by = an 
iterative least-square method to the data points found. The 
data points nowhere devfated from the fitted sphere by more 
than 10 mm, and by tess than 5 mm except near the very edge 
of the image. Such accuracy wilt not usually be obtained 


because of non-uniformity in the paint, shortcomings of the 
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sensing device etc.. For many purposes however less accuracy 
is quite acceptable and for object recognition In particular 
a more important criterion fs most probably that similar 


objects are distorted in similar ways. 


4.1.4 PROBLEMS WITH THE SEQUENT IAL APPROACH: 


It soon become apparent that solving the characteristics 
sequentially had many disadvantages tn the general case, even 
though it works well for lunar topography. The first reason 
is that as the characteristics spread out from the singular 
point, they begin to separate and leave large portions of the 
image unexplored, obtaining only a very uneven sampling of 
the surface of the object (This Is no problem for lunar 
topography since here the solution is not started from the 
singular potnt, but at a place where the spread of the 


characteristics is small), 


With a more parallel approach new characteristics can be 
interpolated as we go along (and some deleted as they 


approach too closely). 


Next we find that the base characteristics (projections of 


the characteristics onto the image) may sometimes cross, 
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Figure 37: Comparison of spread of characteristics for 
typical solutions in case of Tunar topography 
and general case. 


This is not possible if the solution was exact, since it 
indicates that the surface is double-valued or at least that 
its gradient is double-valued, Characteristics may converge 
or diverge froma (singular) point however, Crossing of 
characteristics is really symptomatic of another problem 
which was touched upon when proving the equivalence of the 
five 0.D.E.°s to the P.D.F. : The differential equations for 
p and q must continue to give consistent results with the 
surface calculated - this does happen if the solutton is 
exact, but cannot be hoped for with the noisy data obtained 
from the image. What one would lttke to do Is continuously 
monitor whether the current p andq match with the’ slopes 
obtained by first differences from points on the current and 
neighboring characteristics, This is not possible if the 


characteristics are solved separately and are spreading apart 
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as well. A later section (4.2.2.1) witt explain a method 
used to continously adjust p and qto remain consistent 
(derived from the method explained earlier for finding p and 


q on the initial curve). 


At a very minimum, to avoid embarassment one would like to 
detect when two characteristics approach one another and stop 
one before they cross, This is easy if the solutions are 
carried along in parallel, but involves lengthy comparison 


tests otherwise. 


4.2 A PROGRAM SOLVING THE CHARACTERISTICS IN PARALLEL: 


Once it had been demonstrated that the equations were correct 
and anumerical solution possible fit was decided to write a 
second program which would sample the surface of the object 
more evenly by Interpolating new characteristics when needed, 
Less attention was paid to accuracy in’ the solution while 
attempting to be less sensitive to various noise-effects, At 
the same time an effort to, find a more convenient coordinate 
system produced the much shorter notation and resultant 
equations described in chapter 2, The solution is achieved 
by taking all characteristics one step forward at the sare 


time. 
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surface-section described by Coon [10]. Usually the 
increased complexity tmposed by such an approach can be side~ 
stepped by rather using a smaller step-stze to obtain a finer 


ertd, 


It should be noted that the user of constant size steps along 
the characteristics may produce difficulties on complex 
objects. For even with smooth surfaces the curves. of 
constant arc-distance from the singular potnt (the rings) may 
have cusps, This invalidates the use of difference methods 
on points along these curves (s.a. are used in’ subsection 
4.2.2.1] and 4.2.2.3), No difficulty was experienced with 
images of the objects we experimented with. An alternative, 
which would circumvent this problem, would be the use of 
steps traversing a constant Increment in intensity. This 


would turn the rings Into contours of constant intensity. 


4.2.2 EXTRA PROCESSING POSSIBLE: 


4.2.2.1 SHARPENING ~ UPDATING P AND q: 


We have already described how one can obtain p(t) and a(t) on 


the initial curve by solving the set of non-linear equations 


(subsection 2,3.7): 
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p(t) x(t) + q(t) yy(t) - Ey(t) = 0 


ACr) 4(I,E,G) - b(r?) = 0 


In the proof that the solution of the five ordinary 
different ial equations is also a solution of the original 
partial differential equation, it was stated that the two 
equations for p andq do in fact continue to give the 
derivatives of zwer.t. % and y. When solving a difference 
equation approximation from noisy data we can expect the 
solution for p and q to becore progressively more inaccurate. 
Yet the above pair of equations must hold on any path along 
the surface of the object. In particular one can use them on 


the curve defined by one ring to determine values of p and aq. 


For the initial curve we had the additional difficulty that 
the two equations might have more than one solution and we 
selected one on the basis of some external knowledge (e.g. 
that the object is convex near the singular point). We have 
already assumed that the object is smooth and therefore we 
will have fairly good values for p and q and cannot get into 
this difficulty at non-singular points. Even a simple 
Newton~Raphsen method will suffice to get us more accurate 


values of p anda, 
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Let CARGO?) = By eS og ee ee 
se ecm cp ae bCI,E,G) - b/A 


And suppose: g(ptho, atSq) = h(p+$p, atfg) = 0 


Then tsnoring other than first-order terms we have: 


fp “8 e & (p,q) 


hp by]/\Sa nh (p,q) 


J 


That is: 


Xt Ye $1) #(p,q) 

p Q $a h(p,q) 
Here x, and ye have to be estimated from difference 
approximations. One may not want to apply the full 
correction (Sp, $a). Kore than one iteration will not be 


required since p and q are very close to the correct values, 


4.2.2.2 INTERPOLATION AND CROSSING TESTS: 


When the separation between two neighboring points in a ring 
becomes greater than 1.5 times the step size along the 
characteristic, a new characteristic is interpolated, Its 
X, YY, F, P and aq values are set to the average of its 


neighbors white the backward pointer is set to -l. A more 
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complicated interpolation method can also be used which 
constructs the line of intersection of the tangent planes at 
the two neighboring points. It then finds the point on this 
line closest to the two neighbors and finally uses a point 
half-way between the point determined previously by the 
simoler method and this new point (This, for small angles 
between the tangent planes, is accurate for a_ spherical 
surface), This method does not however add significantly to 


the accuracy of the solution. 


If two neighboring points in a section of a ring come closer 
than 0.7 times the step-size, one is deleted (It is important 
that this factor be less than 0.75, that is, one half of the 
factor used in the interpolation decision, or succesive rings 
on a flat region will have points interpolated on one step, 
onty to be removed on the next, with consequent loss. of 


accuracy). 


Finally one wants to stop neighboring characteristics from 
crossing over each other. Consider the two points a and b on 
one ring and their successors c and don the next. The test 
consists of checking whether c is to the teft of the directed 
Line through bd and whether d is to the right of the directed 
Line through ac (Both tests are needed). If either fails, 


the corresponding characteristic is terminated, causing a 
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Figure 38: The four points used in the crossing test. 


Figure 39: The five neighbors used in determining the 
intensity gradient at P. 


Figure 40: Covering the image with the rasters of points 
read for each solution point. 
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hreak to appear in the ring at that point. The test is 
equivalent to checking whether the Line segment cd falls in 
front of the Line segment ab (and does not cross it). This 
test is applied across short breaks In rings as well, to stop 
neighboring section of the ring from crossing over. each 


other. 


Care has to be taken if the sections of a ring left all fall 
on one side of the singular point, since the break then 
actually encompasses an arc of more than TT and crossing tests 
applied across it will invariable terminate more 
characteristics on either side of it, This can be avoided if 
the crossing test Is not applied to points whose images fall 
too far apart (in terms of the projection of the current 


step-size). 


4.2.2.3 OBTAINING GOOD INTENSITY GRADIENTS: 


To be more noise-jmmune than the previous program, a better 
way had to be found to obtain intensity gradients. Rather 
than use the intensities at asmall raster of points to 
estimate the local gradient, it was decided to use a 
difference approx imation from intensities measured at 


neighboring points, Using as many as possible of the 
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intensities of the point Itself and its five immediate 
neighbors, we can apply a simple least-squares method to 
estimate the gradient, Some of the points may not exist as 
explained previously and the characteristic is terminated If 
less than three points are available or only three which are 
nearly collinear. Suppose the coordinates of the points are 
(xKe Yu) (image coordinate system) and the intensities are by. 
We wish to find by , by: and by:, to minimize the following 


expression: 


Deby xh by YL # by > bq) 
K 


This happens when: 


at oe e o 
SX SX SX by \ (Sb AXE 

o 80 ar ? e 
SX SY => Yk bys J] Sb, Vy 
SxS =! bo Sb, 


From bys and by: we can find by, by and b= by using the 


camera projection equations of an earlter section (2.7). 


For good noise~immunity and some ability to detect surface 
detail indicating that the solution is invalid, the intensity 
for each solution point Is not read from only one image 


point. Small tilted rectangular rasters of points are 


Pare 75] 


established around each point of the solution. The ons axis 
of the rectanste is parallel to the base characteristic at 
that point, and thre sire is adiusted to correspond to. the 
projection on the fimree of a square on the object of side- 
lensth equal to the step-size, The intensity recorded for a 
solution point is the averare of the intensities read for the 
points in this raster and the rem.s./averase is used to make 
the edge-crossing decision. The rasters of all the points in 
the deta-structure almost, but not quite, touch and = taken 
together almost cover the total area of the imase explored, 
This insures that the data is not much affected by pin-holes 
in the photo-cathode of the itmare-dissector and that edre 


crossing can easily he detected, without reducing the 


resolution, 


Both this program and the one discussed jin section 4.1] spend 
more than half their time accessing the imare-dissector. 
Between 20 and 100 intensities are read for each point in the 
solution, and each access takes about ,? to 7.0 milli- 
seconds. A complete solution requires from 1 to 5 minutes 


of real time. 
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4.2.3 A DOZEN REASONS TO TERMINATE A CHARACTERISTIC: 


This is a good place to summarize the reasons for terminating 
the characteristics, The values printed near the end of a 
characteristic (derived from the negative intensity code 


discussed earlier) can be used to index this table. 


1. The characteristic has moved out of the field of view 


of the image-dissector,. 


2. The r.m.s./average for the intensities read in the 
raster has become too great, Indicating overlap of two 
objects or an angular joint on one object or some 


surface detail that Is being missed. 


3, The intensity has becone too low, tndicating a shadow 


region. 


4, X is too large, indicating approach to either another 


singular point or an ambiguity edge. 


5. There are too few neighbors to construct a_e good 


estimate of the local intensity gardient., 
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6. & is too small. QO is the Jacobian of the image 
transformation from zy and zy to 2,’ and zy’ This 
transform becomes singular when = 0, In most cases E 


will become too small before this happens, 


7. A new point was interpolated but both Its nelfghbors 


were terminated before it could get anywhere. 


8. I too small - indicating approach to a shadow edge, 
9. E too small - indicating approach to an edge of the 
object. 


10. This characteristic crossed over a neighboring one, 


11. It was discovered that this point has a backward 
pointer to a nonactive point. This ts really an error 


condition and shouldn’t normally happen. 


12. The intensity is equal to or greater than that measured 
at the singular point, Indicating another singular 


point or ambiguity edge. 


Note that several of these conditions are redundant to ensure 


that even with an inexact solution at least one will fail at 
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the right place. 


4.2.4 OPERATION OF THE PROGRAM: 


4.2.4.7 THE INTEGRATION PROCESS 


First the procram needs to be given such parameters as the 
position of the Light-source, the distance to the object, 
focal length of the lens and the step-size to be used in the 
integration. It then proceeds to find a point of maximum 
intensity (for somo reflectivity functions one needs to 
search for a minimum), This search can be directed to allow 
a choice of one cf several possible maxima. The program then 
assumes that this point of maximum intensity Is a singular 
point and that the object. is convex at this point Cin some 
cases we would Like to assume jit to be concave), After 
constructing an initial curve (a small cirele) around the 
singular point, it preceecs to read the intensities at the 
corresponding imare points, The non-linear equations for »p 


and q on this curve are then solved iteratively, 


All intensities are normalized wyr.t. the intensity at the 
singular point unless the surface has a specular component, 


In the latter case, the intensities on the Initial curve are 
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FIGURE 41A 


FIGURE 41B 


FIGURE 41C 
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Figures 42 A, B, C: Stereo-pairs of solutions produced by 


new SHADE for disc-shaped, spherical and bullet-shaped 
objects (actually spheres with make-up applied). 
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Stereo-pairs of same solutions as in 


Figures 43 A, B, C: 


previous figures, rotated 90°, 


used to establish a normalization value (The specular 
reflectivity is too veritable for use in normalization), It 
is assumed that the initial curve has been chosen targe 
enouph fo fall outside the region ef strong snecutar 


reflection. 


For each step in the parareter s, the following procedure is 


then carried out: 


1. For each point calculate the normal (n), the incident 
vector (r;) and the emittance vector (r,). Fror these 


obtain the derivatives Ty, £y and G,. 


2. Calculate by a, bv, and hence 4 


n? 
~ 


3. Then obtain Foy F, and). 


4, Add (Sx, Sy, fz) Lo (x,y,z) to get the potnt on the next 


tt 


ring for each characteristic. Here (Sx, sy, $z) 


NCE Se Fay pF. +aF, ). 


5. Interpolate new points where the points in the new ring 
are too far apart and delete points where they are tao 
close together. Produce breaks where characteristics 


have crossed over adjacent characteristics, 


12, 


13. 


14, 
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Now read the intensities for all the points. Terminate 
those characteristics with points of very low intensity 


or high r.m.s/average. 


Calculate byt, by: for all those points for which enough 
neighbors exist. From these values obtain bx, by and bz 


by the projection equations, 


Now use n, r- andr, to calculate I,, EF, and G,-. 
~ ~ & ne z : 


~ ~ 


Next use >. b. and >. to calculate b 


Then obtain Fe ee and Fl. 


Add (Sp, $a) to (p,q) to obtain p and q_ for the 
uninterpolated points on the new ring. Here (Sp, $q) = 


AG GR pF) 7 CR ab) sy 
Interpolate p and q for the new points, 


Sharpen up the values for p and q on all points in the 


new ring. 


Garbage-collect various ftems, such as serfes of points 


with negative intensity. 
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15, Stop if no points with positive intensity remain. 


It should be apparent where the vartous tests for terminating 
the characteristics fit into the above schema, The simple 
Euler method for solving the differential equations could be 
replaced by a Runge-Kutta method with increases in running 
time of a factor of two, but little improvement in accuracy. 
The sharpening method, on the other hand, is very cheap and 


contributes substantially to accuracy. 


4.2.4.2 OTHER PROCESSING AVAILABLE 


As explained before, the data-structure is displayed as it is 
generated and can also be viewed from different angles when 
completed, In addition a mode exists where the mapping from 
three-space to the display surface is not performed by the 
projection explained earlier, but a simple map from a 
rectangular area on the image-dissector to a rectangular area 
on the display surface. This is particularly valuable for 
overlaying the solution on an Intensity modulated display of 
what appears in the image. This aids greatly in debugging 
since it is easy to pinpoint such problems as starting the 


solution from an inappropriate maximum in Intensity. 
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A number of other displays can be oroduced to aid in setting 
up the image-dissector. Prodiszous amounts of detailed print- 


out can be generated during a solution process anda more 


parsimonious Listing of the final data is available. It is 
possible to substitute synthetic data (with selectable 
amounts of noise) for the image-dissector input as a 


repeatable way of checking out the program and to tide over 
those days when the image-dissector is being repaired! The 


data can be written to and read from tke disk and tape, 


The stereoscopic display has to be viewed with an appropriate 
pair of lenses which are not always handy. For this reascen a 
routine was provided which produces a contour map from the 
data. This map is produced by first listing the 
intersections of all the Lines in the data structure (from 
point to point in each characteristic, as well as from point 
to point in each ring) with the selected contour planes. The 
intersections are then sorted on contour plane. Within each 
contour plane the foltltowine process is applied repeatedly 


until no points are left: 


Pick a point and find the closest neighbor within a 
“reasonable” distance. “*Reasonable” distance is defined 
to be 1.5 times the step-stze used tn the solution, Now 


another point fs selected closest to the new point also 
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within a reasonable distance and so on until no more can 
be found. The point chosen at each step may not be the 
first in the chain so constructed (which would close the 
loop) unless no other points are available. Atso the 
Line-segments connecting sucessive points may not make 
angles of more than 71/2 with one another. The points 
are removed from the data as they are used in generat ing 
the contour except the very first point (to allow for 


the eventuality of closing the contour), 


The distances are usually weighted with the dot-product of 
the new segment with the previous segment, to give preference 
to continuation of contours in the direction of the last 
segment used, The method generates good contours where’ the 


data is complete and smooth, and does fairly well otherwise. 


4.2.5 INSENSITIVITY TO IMPERFECTIONS IN THE SENSOR: 


This program is not quite as accurate as the one that. solves 
the characteristics serially (mostly because of the simple 
method for solving the differential equations numerically), 
but vastly superior in its behavior when’ faced with noisy 
data. Most of the fmprovement is due to the better way of 


obtaining intensity gradients and to the use of the lateral 
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FIGURE 45A 


FIGURE 45B 
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Figures 47 and 48: Stereo-pairs of solutions obtained for 
the plaster object and the cube with rounded corners. 
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Figure 49: Stereo-pair of side-view of solution obtained 
for the plaster object. 
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Figure 50 and 51: Contour maps of the solutions obtained 
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for the plaster object and the cube with rounded corners. 
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connection between the characteristics. The difference 
approximation for the intensity gradients uses a support area 
about six times as targe as the one used by the least squares 


approximation of the first program, 


Distortions in the imag ing device ‘merely’ produce 
distortions in x and= ey, while non-uniformities in the 
sensitivity will affect p and q and hence 2. The only effect 
of low resolution will be that some edges will not be noticed 


and the solution erroneously continued across them. 


4.3 A NOSE-RECOGNITION PROGRAM: 


To illustrate one use of the shape-from-shading method, it 
was applied to a simple recognition task. Although there is 
great interest in face-recognition (12] (partly because there 
is a practical use for it), it was decided to tackle a sub- 
problem - that of nose~recognition. In principle, face- 
recognition could be carried out by repeating the process 
explained here for not only the nose, but the chin, forehead 
and the two cheeks, Transparencies of noses, rather’ than 
real noses were used because they are always ready and do not 
move during the minute or so it takes to determine the shape. 


To avoid having to determine the reflectivity function for 
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skin as a function of all three angles, special tLighting 
conditions were employed. The light-source was placed near 
the camera and the reflectivity function as a function of the 
incident angle determined from the transparencies taken. 
This meant that no separate determination of the non- 


linearities in the photographic process was needed. 


4.3.1 MODIFICATIONS TO THE BASIC PROGRAM REQUIRED: 


A few minor changes and additions had to be installed in the 
main program for this task. Most prominent amongst these is 
the procedure used to normalize the intensities read from the 
image. Because of the strong specular component of highly 
variable nature, the singular point could not he used for 
this normalization. The specular component In the 
transparencies not only varies from person to person and time 
to time but depends on the exposure used, since it usually fs 
bright enough to saturate the film, Normalization was thus 
carried out w.r.t. an intensity derived from that measured on 
the initial. curve, which was assumed to be outside the 


specular region. 
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Figure 52: Solution obtained for a nose. 
Note gaps left by the breaks caused by 
the nostrils. 
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Figure 54: Four views of solution obtained for a nose. 
(With some hidden lines eliminated). 
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Figure 55: Stereo-pair of solution obtained for a nose. 
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4.3.2 NORMALIZATION PROCEDURE: 


In order to simplify photographing the subjects, it is 
necessary to make some decisions about which factors one is 
going hold fixed and which are to be taken care of by some 
normalization in the program. Altough it is possible to hold 
the head in a standard position by means of a bite-bar, it is 
inconvenient and it ts preferable to let the program take 
care of small head-rotations. The distance from the camera 
to the subject on the other hand is. very easy to determine 
and therefore no normalization of size was used. For 
pictures of the whole head such size normalization would be 
fairly accurate, whereas it cannot be for images of the nace 
alone which does not present sharp features to take 


measurements of. 


The rotational normalization procedure to be described can 
handle quite large (<1TT/6) rotations In both pitch (rotation 
about an ear to ear axis) and roll (rotation about a tip-of- 
nose to back-of-head axis). Yaw (rotation about a_ top-of- 
head to throat axis) is restricted by the requirement that 
almost all of the surface of the nose should be visible. For 
some noses this restricts the rotation to fairly small angles 
- of course this presents no problem when taking the 


photograph. 
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Figure 56: Illustration of rotational normalization 
procedure. 
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Figure 57: Illustration of parameters abstracted from one 
horizontal contour through the nose. 
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Independence of rotation is achieved by means of a routine 
which establishes the orientation of the shape calculated and 
then rotates it into a standard position. In addition the 
parameters in the final comparison procedure where chosen to 
be independent of small remaining errors in the orientations, 
The orientation of the shape calculated is estimated from two 
horizontal contours through the nose, one passing through the 
tip of the nose, the other higher up on the ridge. These 
contours of course are only defined as sequences of points 


where the characteristics and rings pass through each plane. 


The most forward points defined by these contours are 
calculated by fitting a parabola to the three points with 
lowest z coordinates. For each of the two contours we get 
one such forward point, connecting them we obtain a_ Line 
which runs approximately along the ridge of the nose. This 
Line is rotated into a standard position (Lying in the y-z 


plane and teaning TT/6 from the vertical). 


The lower contour (through the tip of the nose) is also used 
to estimate rotation about the vertical axis, The two points 
on this contour at a given distance from the most forward 
point define two angles w.r.t. the z-axis. The desired 


rotation is one half of the difference of these two angles. 
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The three angles so determined are small and can thus” be 
treated independently, The rotation of the shape is 
performed about the center of the spherical cap used to 
determine the initial curve, I.e. a point just inside the tip 
of the nose. The whole process is repeated iteratively three 
times. The errors remaining are almost always less than 0.01 
radian (0.5°), It was found that ustng only the few points 
indicated to determine the rotation was quite satisfactory, 
although better accuracy ts no doubt obtainable If the 


calculation employed averages over several points. 


4.3.3 COMPARISON PROCEDURE: 


After the data has been brought into a standard orientation, 
we would Like to abstract a small number of parameters which 
contain most of the information for comparison purposes. A 
rather arbitrary decision was taken to use estimates of the 
distance of the ridge of the nose from the standard tine (In 
the y-z plane and leaning tr/é6 from the vertical), the width 
of the nose about half-way down to the cheek and the depth of 
the cheek from the ridge of the nose. These quantities where 
measured for each of five horizontal contour planes, the 
lowest through the tip of the nose, the highest a bit below 


the saddle point (the bridge between the eyes). The fifteen 
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Figure 58: The points on the 5 contours used to abstract 
the fifteen values describing this nose. 
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values so obtained are the only data used in the final 


comparison procedure, 


The distance down the side of the nose from where 
measurements of the width of the nose are taken varies with 
the contours, going from some large value for the plane 
passing through the tip of the nose to one-half that value 
for the highest contour, The distance at which the depth of 
the cheek is measured is twice that at which the width of the 
nose is measured and thus also varies from contour to 
contour, The depth of the cheek is the average of the depth 
obtained on the teft side and that obtained on the right, 
The fifteen measurements obtained for each transparency are 
stored in a table together with the number of the 


transparency. 


The purpose of the comparison procedure is to establish if 
any of the stored measurements match those obtained from a 
new transparency, To determine this, a pseudo-distance' is 
calculated (in the 15-dimensional vector space), between each 
stored vector and the new vector. The pseudo-~distance fs a 
weighted r.m.s. of differences fin coordinates [12], where the 
weights are proportional to the standard deviation observed 


for that coordinate. 
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2 
d°= > OfexD / e 
t 


where dis the pseudo-distance, xi and x the components of 
the two vectors, and o; the standard deviation of the ; th 
component. The uncertainty In the depth to the cheek is 
greater than that in the width of the nose, for example, and 
it therefore has a lower weight than the latter. This 


procedure gives a comparison test which is in some” sense 


optimal [172]. 


No doubt other comparison procedures and other choices of 
parameters would have been equally useful; in particular it 
soon become apparent that fewer than 15 parameters would have 
been equally as selective. The point is that once one has 
data as complete as a_ full description of the shape, almost 
any method will work and it is not even necessary to display 


great sophistication in one’s use of statistics. 


4.3.4 RESULTS OF THE NOSE-RECOGNITION PROGRAM: 


15 transparencies of 12 noses were used fin this experiment. 
The pairs of transparencies for the three noses which were 
photographed twice differed in camera to subject distance, 


head rotation and exposure. A total of 30 shapes were 
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calculated, two each for those noses of which only. one 
transparency was available (they differed because of the 
noisy nature of the data). For each shape so determined, the 
rotational normalization was applied and the 15-tuple 
description abstracted, The pseudo-distances between all 


pairs of 15-tuples were then determined. 


The pseudo-distance between 75-tuples averaged to the 


following (the units are about 0.3 mm’s r.m.s.): 


1. Between transparencies of different noses - 10, Crange 
2.4 - 15.3) 

2. Between transparencies of the same nose ~ 2, (range 1.4 
to 2.5) 


3. Between shapes calculated from the same transparency 


1. (range 0.1 to 2.7) 


In all cases the distance from a given shape toa related 
shape was less than a quarter of the distance to a unrelated 
shape. Simply looking for the smallest pseudo-distance (and 
checking whether it is fairly small), thus gives an effective 
recognition procedure for this smatl data-set. It is clear 


that for a much larger data-set unique identification would 
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be more untikely without improving the accuracy’ in tke 
solution and a detailed analysis of which parameters to 
abstract for optimal recognition. It would however always be 
possible to separate out some smalt subset of the total 
stored set of nose-descriptions with very high probability 
that the nose looked for will be in this set. Bledsoe [12] 
uses the ratio of the size of this subset to the size of the 
complete stored set as a measure of the effectiveness of the 


recognition procedure, 


Repeating the operations we described here for the other 
large frontal planes (planes with normal parallel to the z- 
axis), one would obtain a face-recognition procedure. It js 
very likely that the subsets of all stored face-decriptions 
determined by applying the above method to cheeks, chin, 
forehead and nose in turn will have only a small 
intersection. This is not to say that other information 
about the face, not obtainable from the shape-from-shading 
method could not add to the accuracy of such a procedure. It 
must be pointed out that some of the feature points used in 
previous attacks on the face-recognition problem are not 
defined by sharp discontinuities (for example the tip of the 
nose) and could best be obtained from a description of the 


shape. 


Page 1&7 


The restriction about the positioning of the Light-source 
could be removed if one took the trouble to measure the 
reflectivity function in more detail and either recorded the 
positioning of the Light-source or worked out in detail a 
method for finding the single tlight~-source from the shadows 
in the image (which should not be very difficult since we 
know the approximate shape of the object we are looking at), 
The full face-recognition problem was not tackled since it 
would require a great deal more work without further 
illustrating the method of determining shape-from-shading. 
Also it will be noted that the study involved a small set of 
noses - a study with a large data-set would contribute little 


more to the understanding of the method. 


Some of the difficulties encountered when determining the 
shape of noses are perhaps worth mentioning. Firstly, most 
noses are not completely visible from any given point of 
view, Most notably the underside (between the nostrils) is 
frequently not visible, and often a small area on the side of 
the nostrils is also hidden. This forced a_ choice of 
parameters which did not depend on these areas, Naturally 
the information of whether these areas are vistble could in 
itself be useful in the recognition procedure if it could be 


reltably determined, In fact our program does not, because 


of the combination of poor resolution in the image-dissector 
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and the simple-minded edge detector, This could be 
circumvented by placing the Light~source slightly above the 
camera, thus ensuring that there always is a narrow’. shadow 


below the nostrils. 


When the solution is erroneously continued across an edre 
(such as that above the nostrils), a second undesireable 
effect appears because of the sharpening procedure. The 
incorrect coordinates of the points calculated after the edge 
is crossed have some effect on their neighbors due to this 
and thus decrease the accuracy of the solution obtained 


nearby. 


Another problem is that some noses have not. one, but two 
closely spaced tips (probably because of the underlying 
cartilage consisting of two symmetrical parts). This causes 
the characteristic growing from one of these peaks towards 
the other to stop, since it ts approaching another singular 
point. A simple solution consists of chosing the radius of 
the initial curve large enough to completely include both 
singular points. Finnally one finds that some noses 
(particularly those belonging to females) have very low 
ridges near the eyes, making it difficult to determine a 


meaningful value for the width of the nose at that point, 
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It should be noted that the reflectivity function was not 
determined with great precision and no account was taken of 
its variation from person to person, It was not important 
that the shape calculated was very close to that of the nose 
from which the image was taken, but rather that differences 
in the shapes of noses should show up as differences in the 
calculated shapes and that shapes determined from 
transparencies of the same nose should be similar, If the 
images were all produced with the heads in the same 
rotational position, the distortions wouid have made no 
difference at all. For the small head rotations encountered, 
the effect of the relatively minor distortions was very 


small, 


4.4 SUMMARY AND CONCLUSIONS: 


After defining the reflectivity function, an equation was 
found relating the intensity measured in the image ofa 
smooth opaque object to the shape of the object. This 
equation was then shown to be aé_ first-order non-linear 
partial differential equation in two unknowns and _ the 
equivalent set of five ordinary differential equations was 
derived. A number of especially simple cases were 


discussed, in particular applications to lunar topography and 
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the scanning electron microscope. Methods were described for 
obtaining the auxiliary information required (e.g. the 
reflectivity function) and how to avoid the need for an 
initiai Known curve on the object, Of importance too is the 
method demonstrated for continuously updating p and q 


(sharpening) as the solution progresses, 


The half-dozen or so. other depth-cues were ignored here to 
allow a comprehensive treatment of shading. The analytical 
approach to the problem of determining shape from shading was 
developed to demonstrate that an exact solution is possible 
and to determine just what the limitations of this approach 
are, This ts not to say that a more heuristic, approximate 
approach does not have its merits too for certain types of 
objects [14], It was decided to produce a program to allow 
experimentation with the solution method because many ideas 
in the field of artificial intelligence and visual perception 
are of Little value until they can be tried on real data, 
Fortunately an image-dissector was available to provide input 


of image intensities to the computer, 


Two programs were presented, one solving the 0O.D.F.%s for the 
characteristics sequentially, the other in parallel, 
Advantages of the latter approach were found to be several, 


Finally this latter program was adapted to provide input for 
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a nose-recognition procedure. 


It has been made apparent that shadins its valuable as a 
monocular depth-cue altough it may not be as accurate as some 
others, It must be emphasized that no claim is made that 
people employ this depth-cue in the same way. It may be that 
the human visual system does not actually determine the shape 
in three-space and if it does so it ts likely that it uses a 
different method, However there will be many similarities 
between the two systems (e.g. in the errors they make) 


because they utilize the same data. 


4.4.1 SUGGESTIONS FOR FUTURE WORK: 


I}. It would be instructive (but very time consuming) to 
measure many reflectivity functions and see how many 
fall into the pattern of a matt component, approximately 
varying with cos(i), plus a specular component, If it 
could be shown that most real reflectivity functions 
fall into this class, the method presented would be more 
useful stnce it could determine approximate shapes 
without knowing much more about the reflectivity 


function. 
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It may be possible to find more simplifying conditions 
S.a. the ones’ found with certain Lighting conditions, 
positions of the Light-source and special reflectivity 


functions, 


Other solution methods may be found, or modifications to 
the integration method might increase the accuracy. 
Perhaps a difference method on a fixed grid could be 
found which somehow gets around such problems as that of 


ambiguity edges. 


One could study the two related problems of finding the 
reflectivity function, given the shape of the object and 
the Light-source distribution and finding the Cight- 
source distributton given the reflectivity function and 


the shape. 


Further study of certain types of inconsistencies and 
their use is Indicated. Here for example we find the 
problem of deciding whether certain faces in an image of 
several polyhedra could consistently belong to. one 


object. 


Some effort could be directed towards [Implement Ing more 


fully some of the ideas developed theoretically here, 


10. 
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S.a. shadow bridging, handling multiple sources and 


multiple singular points, 


Expanding the nose-recognition program into a full face- 


recognition program would increase its usefulness. 


One could study in more detail how people use the depth- 
cue of shading and how bad animals are at it. Perhaps 
one can get a better clue as to whether people develop a 
three-dimensional model of the object from the shading 
or if they use the shading information in some. other 


Waye 


There are probably a few more loose ends such as the 
problem of how to start the solution if no convex or 
concave singular points are available. Can one do 
anything at all with saddle points (even though they can 
camouflage themselves to be indistinguishable from 


simple convex or concave singular points)? 


In addition to interpolation, is ft  reasonble to 
extrapolate? That is, can one generate new 
characteristics next to a solution sheet to explore new 
areas, In particular when a break appears Ina 


solution surface can it be patched-up later? 


17. 


13. 
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More methods witl have to be found to deal with the 


three-dimensional structure once one has determined it, 


As pointed out earlier, the use of constant size steps 
along the characteristics may not be ideal (remember 
that we can adjust the step-size by choosing a different 
oe One particularly attractive idea would be to use 
steps corresponding to constant Intensity change in the 
image. This would turn the rings tnto constant 
intensity contours, rather than curves of constant arc- 


distance from the singular point. 


Many objects have surfaces whose reflectivity cannot be 
described by a function of three angles, or are so 
specular that our methods are of little avail. One 
might try to discover methods of dealing with such 
objects. Examples are chrome car-bumpers, translucent 


wax, hair and a egtass of water, 
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fail on smooth objects. Until now the image processing of single views has been re- 
stricted to objects which can meaningfully be considered two-dimensional or bounded by 
plane surfaces. 


A number of applications of this method will be discussed including one to Junar 
topography and one to the scanning electron microscope. In both of these cases great 
simplifications occur in the equations. A note on polyhedra follows and a quantitative 
theory of facial make-up is touched upon. 


An implementation of some of these ideas on the PDP-6 computer with its attached image- 
dissector camera at the Artificial Intelligence Laboratory will be described, and also 
a_nose-recognition program. 
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